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Abstract 


The transformation of cold neutral intergalactic hydrogen into a highly ionized warm plasma marks 
the end of the cosmic dark ages and the beginning of the age of galaxies. The details of this pro- 
cess reflect the nature of the early sources of radiation and heat, the statistical characteristics of 
the large-scale structure of the Universe, the thermodynamics and chemistry of cosmic baryons, 
and the histories of star formation and black hole accretion. A number of massive data sets from 
new ground- and space-based instruments and facilities over the next decade are poised to rev- 
olutionize our understanding of primeval galaxies, the reionization photon budget, the physics of 
the intergalactic medium (IGM), and the fine-grained properties of hydrogen gas in the “cos- 
mic web”. In this review we survey the physics and key aspects of reionization-era modeling and 
describe the diverse range of computational techniques and tools currently available in this field. 
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1 Introduction 


About 380,000 years after the Big Bang the tem- 
perature of the Universe dropped for the first time 
below three thousand degrees Kelvin, the primor- 
dial plasma recombined to form neutral hydrogen 
atoms, matter decoupled from the photon field, 
and the radiation presently observed as the cos- 
mic microwave background (CMB) was released. 
The Universe then entered a period devoid of 
any luminous sources cosmologists call the cosmic 
“dark ages”, quietly waiting for gravity to draw 
hydrogen gas into collapsing dark matter halos. 
The formation of the very first stars following 
gas condensation and cooling, some 200 million 
years after the Big Bang, marks the end of the 
dark ages and the beginning of the age of galax- 
ies. These young stellar systems emitted energetic 
UV radiation capable of ionizing local “bubbles” 
of hydrogen gas. As more and more stars lit 
up and the mean density of the diffuse inter- 
galactic medium (IGM) — the main reservoir of 
baryonic material at all epochs — decreased, these 
HII bubbles overlapped and progressively larger 
volumes became ionized and heated. Throughout 
this epoch of reionization and reheating, the all- 
pervading IGM became clumpy under the influ- 
ence of gravity, and acted as both a source of fuel 


for the galaxy formation process, and as a sink for 
absorbing the products of star formation (metals, 
wind mechanical energy, and radiation). 

The first theoretical calculations of inhomoge- 
neous hydrogen reionization (Arons & McCray, 
1970; Arons & Wingert, 1972) assumed a uni- 
form IGM and estimated that quasar ionization 
zones would overlap by redshift z ~ 3. These 
timely studies made immediately clear that: 1) 
initially isolated cosmological HII regions would 
grow in size around the first sources; 2) the thick- 
ness of ionization fronts — of order the mean 
free path of Lyman-continuum (LyC) radiation in 
the neutral IGM, Amfp & 0.27(1 + z)~? proper 
Mpc — would be much smaller than the size of 
these HII regions; 3) HII regions would accumu- 
late faster that they can recombine as the UV 
emissivity of the Universe rose; and 4) at the 
epoch when HII overlap occurred, the photon 
mean free path and the intensity of the ioniz- 
ing background would increase abruptly. Many of 
these salient features are seen in current state- 
of-the-art numerical simulations of this process. 
Today, however, we believe that hydrogen reion- 
ization ended some 0.9 Gyr after the Big Bang 
(Fan et al., 2006) and was likely driven by star- 
forming galaxies rather than quasars (e.g., Shapiro 
& Giroux, 1987; Madau et al., 1999; Gnedin, 2000; 
Robertson et al., 2015). We know that the IGM 
is clumpy and that the optically thick “Lyman- 
limit systems” (LLSs), high density regions that 
trace non-linear and collapsed structures, cause 
the mean free path to LyC radiation to remain 
smaller than the horizon even after overlap (e.g., 
Gnedin & Fan, 2006; Furlanetto & Mesinger, 2009; 
Worseck et al., 2014). And we argue about the 
fraction of ionizing photons that escape from indi- 
vidual galaxies into the IGM and the role played 
by energetic X-ray radiation in altering the ther- 
mal and ionization balance of the early IGM (e.g., 
Madau et al., 2004; Oh, 2001; Venkatesan et al., 
2001; Madau & Fragos, 2017). 

Modern observations of resonant absorption in 
the spectra of distant quasars show that reioniza- 
tion was nearly completed at z ~ 6 (Fan et al., 
2006; McGreer et al., 2015; Becker et al., 2015). A 
very early onset of reionization is disfavoured by 
the final full-mission Planck CMB anisotropy anal- 
ysis, which prefers a late and fast transition from 


a neutral to an ionized Universe, and by a vari- 
ety of astrophysical probes (Schroeder et al., 2013; 
Banados et al., 2018; Ouchi et al., 2010; Schenker 
et al., 2014; Onorbe et al., 2017; Villasenor et al., 
2022). Despite a considerable community effort, 
however, many key aspects of this process, such as 
the very nature of the first sources of UV radia- 
tion, how they interacted with their environment, 
the back-reaction of reionization on infant galax- 
ies and its imprint on their luminosity function 
and gas content, the impact of patchy reioniza- 
tion on observables, and the thermodynamics of 
primordial baryonic gas, all remain highly uncer- 
tain. Even a complete knowledge of the sources 
and sinks of UV and X-ray light in the early Uni- 
verse, which astrophysicists by no means possess, 
would not automatically translate into a detailed 
understanding of the reionization process as a 
whole, as this ultimately requires extremely chal- 
lenging cosmological numerical simulations that 
self-consistently couple all the relevant physical 
processes — dark matter dynamics, gas dynam- 
ics, self-gravity, star formation/feedback, radiative 
transfer, non-equilibrium ionizations and recombi- 
nations, chemical enrichment, heating and cooling 
— over a huge range of scales (e.g., Gnedin, 2014; 
Ocvirk et al., 2016; Pawlik et al., 2017; Finlator 
et al., 2018; Kannan et al., 2022). Over the next 
decade, the large wavelength coverage, unique sen- 
sitivity, and manifold spectroscopic and imaging 
capabilities of the James Webb Space Telescope 
(JWST) will allow the discovery of star-forming 
galaxies out to redshift 15, the tracing of the 
evolution of the cosmic star-formation rate den- 
sity to the earliest times, the determination of 
the cross-correlation of galaxies with the opacity 
of the early IGM, the study of the highly ion- 
ized near-zones around z > 6 quasars, and will 
provide tight constraints on both the reionization 
history and the contribution of different sources to 
the ionizing photon budget of the Universe. These 
observations will be complemented by a number 
of massive datasets from several existing facilities, 
such as the Dark Energy Spectroscopic Instrument 
(DESI), the Subaru Hyper Suprime-Cam (HSC), 
the Extremely Large Telescope (ELT), the Ata- 
cama Large Millimeter Array (ALMA), and the 
Hydrogen Epoch of Reionization Array (HERA). 
Even more advanced observational capabilities 
will become available in the very near future, such 
as 30-meter class telescopes, the Euclid satellite, 


the Roman Space Telescope, intensity mapping 
instruments such as SPHEREx, CCat-p, TIME, 
CONCERTO, or COMAP, and, in a more distant 
future, the ultimate 21 cm machine, the Square 
Kilometre Array (SKA). It is generally believed 
that the coming years will mark the beginning 
of a golden age of early Universe and epoch of 
reionization science. 

It seems then timely to survey in a “liv- 
ing” document the key aspects of reionization-era 
modeling and the diverse range of computational 
techniques and tools currently available in this 
field, with the aim of providing a comprehensive 
road map for guiding and interpreting forthcom- 
ing observations of the early Universe. A note 
on terminology. In the remainder of this review, 
we call analytical models of reionization those 
based on a set of equations for some statistical 
description of the process, such as the averaged 
neutral fraction, even in the absence of a closed 
form mathematical solution. By contrast, numeri- 
cal methods of retonization are those that provide 
an actual realization of the hydrogen density field 
in a 3-dimensional volume, starting from initial 
Gaussian density fluctuations. A full hydrodynam- 
ical simulation is thus a numerical model, but so 
is the “semi-numerical” realization of the HII field 
provided by, say, the 21cmFast code (Mesinger 
et al., 2011). 

In the spirit of a “living” article that is kept 
up-to-date, we have tried to structure this review 
in a way that maximizes the visibility of the 
latest findings for the reader, critically evaluate 
existing work, place it in a meaningful context, 
and suggest areas where more research and new 
results are needed. Given the limited space avail- 
able, it is impossible to provide here a thorough 
survey of the huge community effort behind reion- 
ization and IGM studies without leaving out 
significant contributions or whole subfields. For 
additional material, we therefore refer the reader 
to a number of recent complementary reviews 
on these topics (Meiksin, 2009; Choudhury, 2009; 
Trac & Gnedin, 2011; Zaroubi, 2013; Lidz, 2016; 
McQuinn, 2016; Dayal & Ferrara, 2018; Robert- 
son, 2021). Unless otherwise stated, all results 
presented below assume a flat cosmology with 
parameters (Qy7,Qa,2,) = (0.3,0.7,0.045), a 
Hubble constant of Hp = 70 km s~t Mpc™!, and 
a primordial baryonic gas of total hydrogen and 
helium mass fractions X = 0.75 and Y = 0.25. We 
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Fig. 1 Cosmic Reionization Fundamentals. After recombination at z ~ 1100, hydrogen in the IGM remained neutral until the first 
stars and galaxies formed at z ~ 15. These primordial systems released energetic UV photons capable of ionizing local bubbles of 
hydrogen gas. As the abundance of these early galaxies increased, HII bubbles overlapped and progressively larger volumes became 
ionized. The hydrogen reionization process completed at z ~ 6, approximately 1 Gyr after the Big Bang. Observations suggest that 
the double reionization of helium occurred at a later time (e.g., Kriss et al., 2001; Shull et al., 2004; Zheng et al., 2004; Shull et al., 
2010; Worseck et al., 2016), when hard UV-emitting quasars and other active galactic nuclei (AGNs) became sufficiently abundant. 


(From Robertson et al. (2010); used by permission.) 


shall also use spectroscopic notation to indicate 
the ionization states - HI and HII for neutral and 
ionized hydrogen, and Hel, Hell, and HelIII for 
neutral, singly ionized, and doubly ionized helium, 
respectively. 


2 The Physics of Intergalactic 
Gas 


The IGM is the storehouse of the bulk of cosmic 
baryons at high redshift. It forms a filamentary 
network that hosts the formation of galaxies and 
gives origin to a “forest” of Lyman-a absorption 
lines, the signature of a highly fluctuating medium 
with temperatures characteristic of photoionized 
gas (see, e.g., Meiksin, 2009; McQuinn, 2016, and 
references therein). It is dark matter that pro- 
vides the backbone of large-scale structure in the 
Universe, a web-like pattern present in embry- 
onic form in the overdensity motif of the initial 
fluctuation field and only sharpened by nonlin- 
ear gravitational dynamics (Bond et al., 1991). 
The Lyman-a forest traces this underlying “cos- 
mic web” on scales and at redshifts that cannot 
be probed by any other observable. Because of its 
long cooling time, diffuse gas at z ~ 2 — 5 retains 
some memory of when and how it was reheated 
and reionized at z 2 6 (Miralda-Escudé & Rees, 
1994). The physics that governs the properties of 


the IGM throughout these epochs remains simi- 
lar, as the evolving cosmic UV emissivity and the 
transfer of that radiation through a medium made 
clumpy by gravity determine both the details of 
the reionization process and the thermodynamics 
of the forest. 


2.1 Photoionization of hydrogen 
and helium 


Photons with frequencies v > vyyz will ionize inter- 
galactic hydrogen at a rate per neutral hydrogen 
atom of 


oe) Jy 
Tur = an | hie our(v)dy, (1) 


HI 


where oy1(v) is the bound-free absorption cross- 
section from the ground state and hvyy = 13.6 eV 
is the ionization potential of the hydrogen atom. 
Here, J, is the angle-averaged specific intensity, 


Ress = : dQ1(t,.x,n), (2) 


TT 


where I,(t,x,n) is the intensity of the radiation 
field at cosmic time ¢ and position x, measured 
along the direction n. Protons will radiatively 
capture free electrons at the rate per proton 


neon (T), where 
ese er ame 4a 


is the Case-A temperature-dependent radiative 
recombination coefficient onto all atomic levels 
of hydrogen, T, is the gas temperature in units 
of 104K, and ne is the proper electron density. 
There exist more accurate fits for aj{;, but this 
approximate form typically suffices for analytical 
calculations. Case-A recombination is the correct 
choice for modeling reionization, as H-ionizing 
photons emitted in radiative transitions from the 
continuum directly to the 1s level are typically 
redshifted below threshold before being absorbed 
by the IGM (Kaurov & Gnedin, 2014). In terms 
of the ionization fractions yy = nyi/ny and 
Luu = Nut /ny, the rate equation for low density 
photoionized gas can be then written as 


daxyy 
dt 


A 
= —ryi ly + NeXHNAYT, (4) 


and similarly for the ionization states of helium 
THel = NHel/NHe, THell = NHell/NHe, and Lyerl = 
MHeltt/NHe: 


cil = -ayel er + netuen (After (5a) 
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where d/dt is the Lagrangian derivative and the 
recombination coefficient of helium includes a 
dielectronic capture term aff.;; from Hell into 
Hel (Aldrovandi & Pequignot, 1973). The rate 
equations above neglect collisional ionizations by 
thermal electrons and must be supplemented by 
three closure conditions for the conservation of 
charge and of the total abundances of hydrogen 
and helium. 


2.2 Heating and cooling 


The temperature of an unshocked gas element of 
overdensity A = p,/(pp), where pp is the gas mass 
density, exposed to an ionizing UV photon field, 


redshift 


Fig. 2 Temperature evolution (Tp in units of 10* K, red 
curve) of an IGM gas parcel at the mean cosmic density 
(A = 1). In this example, a constant, power-law (a = 1) 
UV background flux, normalized to yield Py; = 1071? s~ 
and not energetic enough to photoionize HelI, is ’switched 
on’ at 6 < z < 6.5. Also depicted are the hydrogen pho- 
toheating (Hur, dash-dotted blue curve), recombination 
(Arad; dashed brown), Compton (Ac, dashed turquiose), 
and expansion cooling (Kexp, dashed purple) rates per 
H-atom (i.e. H; = nj, H;/ny and A; = A;/nzz), in units of 


10-275 erg sl. 


obeys the equation (e.g. Hui & Gnedin, 1997) 


dT _2TdA , T dp 
a an ae yu dt (6) 
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where H is the expansion rate, yu is the mean 
molecular weight, n is the proper number den- 
sity of all species including electrons, A is the 
total cooling rate per unit physical volume, and all 
other symbols have their usual meaning. The total 
photoheating rate H is summed over the species 
i =HI, Hel, and Hell, H = 50, n; Hi, where 


A; = an iy — hy;) ojdv. (7) 
vy, pw 


In equation (6), the first two terms on the right- 
hand side account for adiabatic cooling, the third 
describes the change in the internal energy per 
particle resulting from varying the total number 
of particles, and the fourth is the net heat gain 
or loss per unit volume from radiation processes. 
Suitable tabulations of the cooling rates can be 
found in Hui & Gnedin (1997) and Theuns et al. 


(1998). Equation (6) must be integrated in con- 
junction with the rate equations (4-5) describing 
the evolving fractional abundances of the three 
species. Figure 2 shows the evolution of the tem- 
perature Tp for an IGM gas element at the mean 
cosmic density (A = 1) photoionized by a constant 
UV background that is suddenly ‘switched on’ 
between redshifts 6 and 6.5. Also depicted are the 
corresponding hydrogen photoheating, recombina- 
tion, Compton, and expansion cooling rates per 
H-atom. In this example, the power-law (a = 1) 
radiation flux is not energetic enough to produce 
Helll. The gas temperature reaches a peak of 
about 15,000 K as the IGM is flash photoionized 
and the abundances of HI and Hel quickly drop 
by 3 orders of magnitude. After reionization, T’ is 
determined by photoionization heating, adiabatic 
cooling, and to a smaller extent by Compton cool- 
ing. Assuming a power-law background spectrum 
of the form J, = Jo(v/vyr)~°%, and approximat- 
ing the photoionization cross-section as oy; = 
o0(v/va1)~°, we can rewrite the photoionization 
and photoheating rates of hydrogen as 


A4traoJo 
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(8) 
and 
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Here, the last expression is valid in a highly ionized 
gas that satisfies the equilibrium condition 


ny 


nal yt = ; (10) 


trec 
where tree = (n-ajq},)~+ is the hydrogen recombi- 
nation timescale. When this condition is satisfied, 
the photoheating rate does not depend on the 
amplitude of J, but only on its spectral shape, 
which opens up the possibility of constraining the 
nature of ionizing sources from measurements of 
the post-reionization temperature of the IGM. 
More realistically, it is believed that cosmic 
reionization was an extended process that started 
around z ~ 10, and that the cosmic web actu- 
ally experienced two main reheating events. Figure 


Redshift z 


Fig. 3 Redshift evolution of the IGM temperature To from 
the best-fit (black line) and 95% confidence interval (gray 
band) model for the UV background photoheating and pho- 
toionization rates of Villasenor et al. (2022). The figure 
displays the epoch of cosmic expansion cooling at z < 6 fol- 
lowing the reionization of hydrogen by massive stars. The 
onset of helium reionization by AGNs at z ~ 4.5 initiates 
a second epoch of heating that ends at z ~ 3 when Hell 
reionization completes. A second epoch of cooling due to 
cosmic expansion then follows. The data points with error 
bars show the values of To inferred from observations of 
the Lyman-a forest by Bolton et al. (2014); Hiss et al. 
(2018); Walther et al. (2019); Boera et al. (2019); Gaikwad 
et al. (2020, 2021). The ACDM simulations used for sta- 
tistical analysis assume a spatially homogeneous ionizing 
background. (Image courtesy B. Villasenor.) 


3 shows the thermal history of the IGM at 
z < 6 as recently inferred using a suite of more 
than 400 high-resolution cosmological hydrody- 
namics simulations of structure formation that 
self-consistently evolve a wide range of physically- 
motivated ionization and thermal histories of the 
IGM (Villasenor et al., 2022). A statistical com- 
parison with the observed transmitted flux power 
spectrum of the hydrogen forest over the redshift 
range 2.2 < z < 5 and with the Lyman-a opac- 
ity of intergalactic HeII over the redshift range 
2.4 < z < 2.9 produces the best-fit temperature 
history depicted in the figure. A first peak in To, 
due to hydrogen reionization at z = 6, is followed 
by an epoch of cooling due to adiabatic expan- 
sion until the onset of helium reionization from 
extreme UV radiation emitted by AGNs. The ion- 
ization of helium leads to a second increase of 
the temperature until Hell is fully ionized (z ~ 


3), which is then followed by a second period of 
cooling from Hubble expansion. 


2.3 Secondary ionizations by X-rays 


It is generally recognized that energetic X-ray 
photons emitted by active galaxies and/or X-ray 
binaries may have significantly altered the thermal 
and ionization balance of the early neutral IGM. 
Since only highly energetic photons result in more 
than one ionization, this process is only important 
for radiation emitted by AGNs and X-ray binaries, 
while stars hardly emit any photons that are able 
to ionize hydrogen more than once. Unlike stel- 
lar UV photons, X-ray photons can easily escape 
from their host galaxies, quickly establish a homo- 
geneous X-ray background, make the low-density 
IGM warm and weakly ionized prior to the epoch 
of reionization breakthrough, set an entropy floor, 
reduce gas clumping, and promote gas-phase Hz 
formation (e.g., Oh, 2001; Venkatesan et al., 2001; 
Madau et al., 2004). Photoionizations of hydrogen 
and helium by X-rays produce fast photoelec- 
trons that can collisionally excite and ionize atoms 
before their energy is thermalized (e.g., Shull & 
van Steenberg, 1985). To account for the effect of 
X-ray “secondary” ionizations, the standard pho- 
toionization rates of the neutral species HI and 
Hel must be “enhanced” as (see Ricotti et al., 
2002; Kuhlen & Madau, 2005; Madau & Fragos, 
2017) 


es 
lw =P + 2 ey see 11 
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Here, the sum is taken over 1 =HI, Hel, and 
Hell, and H; are the standard photoheating rates 
per atom (eq. 7). The above expressions are valid 
in the high-energy limit, in which the fractions 
of the primary electron energy that go into sec- 
ondary ionizations of HI and Hel, leet and ce 
are independent of the electron energy (Shull & 
van Steenberg, 1985). Each correction term for 
secondary ionizations is proportional to n;, the 
number density of the species 7 that is the source of 
primary photoelectrons. Energy losses from elec- 
tron impact excitation and ionization of HeII are 
small and can be neglected, i.e. one can safely 


assume oan = Tyen. Collisional excitations of 
helium also generate line photons that can ionize 
atomic hydrogen, but this contribution amounts 
to less than 5% of the direct ionizations by sec- 
ondaries. Electrons with energies < 10.2 eV are 
unable to interact with any atoms or ions, so all of 
the energy is deposited as heat. Since a fraction of 
the energy of the primary photoelectron goes into 
secondary ionizations and excitations, all photo- 
heating rates must be reduced accordingly by the 
factor fheat, _ 
Hi; = feat H;. (13) 
Thus, the reduced energy deposition rate from 
species 7 is related to the excess photoionization 
rate of species j by the factor (fheathv;/f2,,)- 
Monte Carlo simulations have been run to 
determine the fractions fH, fel, and fneat as a 
function of electron energy and hydrogen ionized 
fraction xy (Shull & van Steenberg, 1985; Valdés 
& Ferrara, 2008; Furlanetto & Stoever, 2010).1 
Here, we use the fitting functions provided by 
Shull & van Steenberg (1985), which apply in the 
limit EF > 0.3 keV, 
HI _ 9. 3908(1 — a9,4092)1-7592 
ion, = 0.0554(1 — afr '*)*°°, = (14) 
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For photon energies above 0.2keV, the ratio of 
the HI/Hel photoionization cross sections (see 
Section 2.6) drops below 4 per cent: since the pri- 
mordial ratio of helium to hydrogen is about 8 
per cent, the photoionization and heating rates are 
dominated by helium absorption, which exceeds 
the hydrogen contribution by ~ 2: 1. While Hel 
is the main source of hot primary photoelectrons, 
however, it is HI that is more abundant and there- 
fore undergoes the bulk of secondary ionizations. 
A primary nonthermal photoelectron of energy 
FE = 0.6 keV in a medium with ionization frac- 
tion xy = 0.01 creates about a dozen secondary 
ionizations of hydrogen (corresponding to a frac- 
tion fH! = 0.29 of its initial energy), and only one 
secondary ionization of Hel (fie! = 0.045), there- 
fore depositing in the process 34% of its initial 
energy as secondary ionizations, 37% as heat via 


!The energy deposition fractions are insensitive to the total 
hydrogen density, which enters only through the Coulomb log- 
arithm in the electron scattering cross-section (Shull & van 
Steenberg, 1985). 


Coulomb collisions with thermal electrons, and the 
rest as secondary excitations. Once the IGM ion- 
ized fraction increases to xy 2 0.1, the number 
of secondary ionizations drops (fH#1+ fH#e! < 0.19), 


and the bulk of the primary’s energy goes into heat 
(fneat 2 0.64). 


2.4 Radiative transfer in an 
expanding universe 


Ignoring scatterings into and out of the beam, 
the evolution of the specific intensity I, (t,x,m) in 
an expanding Universe is described by the radia- 
tive transfer equation (e.g., Peebles, 1971; Gnedin 
& Ostriker, 1997; Abel et al., 1999; Petkova & 
Springel, 2009): 


101, n dl, H(t) ( dl, 
c Ot | a Ox c (v Ov 31.) (15) 
=—kyl, + Sy, 


where x is the comoving spatial coordinate, n is a 
unit vector along the direction of propagation of 
the ray, S, is the source term or emission coeffi- 
cent (with units erg s~' cm~? Hz~! sr~*), and k, 
is the absorption coefficient — the product of an 
atomic absorption cross-section and the number 
density of atoms that can interact with photons 
of frequency v. Equation (15) is essentially the 
Boltzmann equation for the phase-space density 
I, /v, with sources and sinks on the RHS. It can 
be recognized as the classical equation of radiative 
transfer in a static medium with two modifica- 
tions: the denominator a(t) = 1/(1 + z) in the 
second term, which accounts for the changes in 
path length along the ray due to cosmic expansion, 
and the third and fourth terms, which account 
for cosmological redshift and dilution of intensity. 
Comparing the third with the second term, one 
recovers the classical transfer equation in the limit 
when the scale of interest Z is much smaller than 
the horizon size, L < c/H(t). Applying the space- 
and direction-averaging operator ( ) to Equation 
(15), we obtain the following equation for the 
mean specific background intensity J, = (J,): 


OJ, OJ, = a Cc 
A a. Kv Vv a2) 
At (v av 3h) Chy J, +7 (16) 
where €, = 47(S,) is the mean proper volume 


emissivity. The above equation admits an integral 


solution 


I) =— | * dtlew(t) [2] eres, (17) 


where v’ = va(t)/a(t’), 
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7(v,t',t) =c i dt" Ryn (t"), (18) 
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and v” = va(t)/a(t”). A packet of photons will 
travel a proper mean free path Amfp(v) = 1/R 
before suffering a 1/e attenuation. 

The mean opacity &, in Equation (16) is 
defined as ky = (kK I,)/J, (Gnedin & Ostriker, 
1997). This is not the space average of the absorp- 
tion coefficient, since it is weighted in the radiative 
transfer equation by the local value of the specific 
intensity I,. The relation &, = («,) holds only 
in the limit of a uniform and isotropic ionizing 
background, or when «, and J, are indepen- 
dent random variables. The former is a reasonable 
approximation during the late stages of reion- 
ization, when ionized bubbles expand into low- 
density regions under the collective influence of 
many UV sources (Iliev et al., 2006a), the pho- 
ton mean free path is almost insensitive to the 
strength of the local radiation field, and varia- 
tions in the LyC opacity are relatively modest. 
The assumption of a nearly uniform UV _ back- 
ground, however, breaks down during the early 
stages of reionization, when rare peaks in the den- 
sity field that contain more absorbers as well as 
more sources ionize first, and the photon mean 
free path is comparable to or smaller than the 
average source separation. A detailed modeling of 
the spatial correlations between sources and sinks 
of ionizing radiation can only be achieved using 
radiative transfer simulations or semi-numerical 
schemes, which we describe in the rest of this 
review. 

When the photon mean free path is much 
smaller than the horizon size, cosmological effects 
such as source evolution and frequency shifts 
can be neglected. In this local source T — oo 
approximation, the mean specific intensity relaxes 
to 

Ey 


Andy ~ —, (19) 


Ky 


and only emitters within the volume defined by 
an absorption length 1/%, contribute to the back- 
ground radiation field (Zuo & Phinney, 1993; 
Madau et al., 1999). While, for a given ionizing 
photon emissivity, the local source approximation 
can overestimate the background intensity close 
to the hydrogen Lyman edge at z ~ 2 (this is 
because a significant fraction of the emitted LyC 
photons at these epochs gets redshifted beyond the 
Lyman limit and does not contribute to the ioniz- 
ing flux), it quickly approaches the exact value of 
J, at redshifts z > 5 (Becker & Bolton, 2013). 


2.5 The source term 


Both because of finite resolution and neglected 
physics, it is customary to separate the transport 
of radiation in the IGM (and in the circumgalac- 
tic medium or CGM) from that in the interstellar 
medium (ISM) of the source population by intro- 
ducing an escape fraction parameter fes.c — the 
fraction of LyC photons that escape the star form- 
ing or nuclear activity regions and leak into the 
IGM. The escape of ionizing radiation from each 
individual galaxy varies with propagation direc- 
tion, time, and scale and depends on the detailed 
spatial distribution of neutral gas, stars, and dust 
in the ISM as well as the clumpiness of halo gas. 
It is often convenient to sweep these complexities 
into a single parameter, which should be thought 
of as a global average over direction and (some 
of) host galaxy properties. If we denote then with 
¢o(L,,t)dL, the source rest-frame UV luminosity 
function (LF), i.e. the comoving space density of 
emitters with monochromatic luminosities in the 
range (L,,L, + dL,), the mean proper ionizing 
emissivity can be written as 


eAt) = al)" [ fesc,»(t)O(Ly,t)dL,. (20) 


The rest-frame UV (1500 A) LF of galaxies is well 
measured down to Myy = —15 mag (Bouwens 
et al., 2022), and has a steep faint-end slope, 
o(L,) «x LS with a ~ —2. A slope exactly equal to 
—2 implies an equal contribution to the radiation 
emissivity per unit logarithmic interval of galaxy 
luminosity. Thus, each unit magnitude interval 
between the exponential cutoff of the Schechter 
function M* + —20 mag and —15 mag contributes 
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Fig. 4 Redshift evolution of the mean emission rate of 
hydrogen ionizing photons into the IGM inferred from 
measurement of the evolving galaxy (assuming a constant 
fesc = 0.2) and quasar (fesc = 1) LFs. To highlights 
uncertainties in the faint-end cut-off of the galaxy LF, the 
three dashed curves show nion for galaxies brighter than 
Myy = —18,-15, and —13 mag. Similarly, the hatched 
region encompasses the range of values obtained by inte- 
grating the LF down to Myy < —18 or Muy < —21 
mag. The shaded regions show the 68 per cent confidence 
regions for Nion inferred using different observational con- 
straints on the reionization timeline: 1) the CMB electron 
scattering optical depth and the Lyman-a and Lyman-{ 
forest dark pixel fraction (gray shade); and 2) same as 1) 
including additional measurements of the hydrogen neu- 
tral fraction from the Lyman-a damping wings observed in 
quasar spectra (purple shade). (From Mason et al. 2019b, 
used by permission. ) 


equally about 20% of the total UV photon out- 
put. Unless such a steep faint-end slope continues 
down to extremely faint magnitudes of —10 mag 
or so, the implication is that we may have already 
observed the bulk of the galaxies responsible for 
the reionization of cosmic hydrogen. A challenge 
to this line of argument may come from the fact 
that rest-frame UV luminosities do not directly 
translate into ionizing luminosities — the two are 
related by the relative escape fraction, 


L; = Sescion RL F 
eee °  -1500A’ 
esc, 1500A 


(21) 
where R is ratio of the intrinsic, non-attenuated 


ionizing and 1500 A luminosities. Here, f aan 
esc, 


depends primarily on the dust abundance and 
spatial distribution, while fescion is thought to 
depend primarily on the distribution of neutral 
gas from the scales of the disk ISM to those char- 
acterizing the structure of molecular clouds (and, 


perhaps, on the abundance of “runaway” massive 
stars that travel through interstellar space with 
anomalously high velocities) (Gnedin et al., 2008; 
Wise & Cen, 2009; Yajima et al., 2009, 2011; Con- 
roy & Kratter, 2012; Alvarez et al., 2012; Mitra 
et al., 2013; Kim et al., 2013a,b; Kimm & Cen, 
2014; Roy et al., 2015; Ma et al., 2015; Trebitsch & 
Blaizot, 2016; Xu et al., 2016; Kimm et al., 2017; 
Howard et al., 2017, 2018; Sumida et al., 2018; 
Katz et al., 2018; Kimm et al., 2019a,b; Ma et al., 
2020a; Yoo et al., 2020; Yeh et al., 2022). Simu- 
lations of radiative feedback in molecular clouds 
support the view that the escape fractions are 
most sensitive to the short temporal and spatial 
scale processes in molecular clouds rather than 
the global properties of the host galaxies (Howard 
et al., 2017, 2018; Kim et al., 2018, 2019, 2021). 
If that conclusion was fundamentally incorrect, 
variations of the escape fractions could impact 
the contribution of different galaxy types to the 
total ionizing photon budget. While the escape 
fraction is a function of frequency or a spectral 
band, in order to simplify notation we shall use 
below a common shorthand fs. to indicate the 
escape fraction of ionizing radiation, i.e. hereafter 
fesc = fesc,ion: 

Galaxy-dominated scenarios that satisfy many 
of the observational constraints can be con- 
structed if one extrapolates the UV LF to suf- 
ficiently faint magnitudes (Muy ~ —13) and 
assumes an escape fraction of fose ~ 10% at z > 6 
(e.g., Madau et al., 1999; Robertson et al., 2010; 
Finkelstein et al., 2012; Haardt & Madau, 2012; 
Shull et al., 2012a; Kuhlen & Faucher-Giguére, 
2012; Robertson et al., 2015; Bouwens et al., 
2015). In AGNs, fese is observed to be between 
44 and 100% (Grazian et al., 2018). There is a 
vast literature on measuring escape fractions from 
star-forming galaxies, and surveying it is beyond 
the scope of this review. In general, however, 
population-average escape fractions from galaxies 
of around 10% may not be unreasonable. Another 
often-used quantity is 


ition (t) = a(t)? / . ol) dv, (22) 


yi 


the mean emission rate into the IGM of ionizing 
photons for species i =HI, Hel, and Hell per 
unit comoving volume — not to be confused with 
the production rate, as the latter does not account 
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for local absorption and hence does not depend on 
the escape fraction. Figure 4 (from Mason et al. 
2019b) shows the uncertain H-ionizing photon 
production rate jo, inferred from measurements 
of the evolving galaxy and quasar LFs. The shad- 
ings mark the 68 per cent confidence regions 
for Mion inferred using a number of observational 
constraints on the reionization timeline. 

Recombination radiation from a clumpy, ion- 
ized IGM can also provide a non-negligible contri- 
bution (up to tens of percent, Haardt & Madau, 
2012) to the ionizing photon field. The main 
processes that contribute to the H-ionizing emis- 
sivity include recombinations from the continuum 
to the ground state of HI, Hel, and Hell, as 
well as HelII Balmer, two-photon, and Lyman- 
a emission (Haardt & Madau, 1996; Sokasian 
et al., 2003; Haardt & Madau, 2012). The mean 
proper emissivity from a generic recombination 
process of species i can be written as e, = 
hyvg(v)(ainen:), where $(v) is the normalized 
emission profile and a; is the relevant recombina- 
tion coefficient. The emission profile of free-bound 
recombination radiation can be computed via the 
Milne detailed-balance relation, which relates the 
velocity-dependent recombination cross section to 
the photoionization cross section, while a delta- 
function profile is sufficient for bound-bound tran- 
sitions. 


2.6 Continuum absorption 


In a homogeneous primordial IGM, the mean free 
path for LyC absorption of a photon of frequency 
v at redshift z is given by 


Amfp(V, Z) = (= ne] . (23) 


where n; and o; are the proper number densities 
and photoionization cross-sections measured at z 
and v of species i = HI, Hel, and Hell. Accurate 
fits to the bound-free absorption cross-sections 
from the ground state for many astrophysically 
important atoms and ions are given by Verner 
et al. (1996). There is also an analytical for- 
mula for hydrogen (and, by extension, for any 
hydrogenic atom) that is exact in the Bohr atom 


approximation (Osterbrock & Ferland, 2006): 


e(4-4 tan! «/a) 


1 = e727 /x Fi (24) 


ont(v) = 90 (Vr/v)* 


where o9 = 6.30 x 107'cm? and x = 
/ (v/v) — 1. The cross-sections drop off rapidly 
with increasing photon energy, a(v) « v~*-° in the 
X-ray regime, for all atoms and ions, so X-ray pho- 
tons can penetrate significantly further into the 
IGM than those with UV energies close to the ion- 
ization thresholds. For crude estimates, the effec- 
tive bound-free cross-section per hydrogen atom 
in a neutral IGM can be approximated as opp = 
ont + fueTHe ~ (3.3 x 10773 cm?) (hy/keV)—3-? 
to an accuracy of better than 20% in the photon 
energy range 0.07 < hy < 8 keV, where fHe = 
Y/(4X) = 0.083 is the helium to hydrogen ratio 
by number. Over this energy range, helium dom- 
inates photoabsorption over hydrogen by a factor 
of 1.4-2.6. The mean free path for photoelectric 
absorption in a neutral IGM is then 


Amfp(¥, 2) © S2Mpe (me) ean (25) 


keV 10 


in physical units. Because of the strong energy 
dependence of the cross-section, photons with 
hv < 1 keV are absorbed closer to the source 
galaxies, giving origin to a fluctuating UV/soft X- 
ray background. More energetic photons instead 
permeate the Universe more uniformly, heating 
the low-density IGM far away from galaxies. The 
pre-reionization Universe is optically thin to all 
radiation for which Amp > c/H, ie. photons 
above ~ (1.6 keV) [(1+.2)/10]°-4” have a low prob- 
abilitiy of absorption across a Hubble volume 
(McQuinn, 2012). 

The IGM, however, is known to be very 
clumpy. Its diffuse, close-to-mean-density compo- 
nent is observed to be highly photoionized at 
z S 6, with only a small residual amount of neu- 
tral hydrogen set by the balance between radia- 
tive recombinations and photoionizations from a 
nearly uniform UV radiation background, and 
provides negligibly small HI photoelectric absorp- 
tion. The continuum optical depth is instead 
dominated by the LLSs, high density regions in 
the outskirts of galaxies (Fumagalli et al., 2011; 
Shen et al., 2013) that occupy a small portion of 
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Fig. 5 Volume-weighted PDF Py(logA) of the 
baryon overdensity A from cosmological hydrodynam- 
ics simulations in which a uniform UV background is 
switched on at z = 9 (Bolton & Becker, 2009). The 
curves shows the PDF at redshifts z = 7 (magenta 
solid line), z = 5 (blue dashed line), and z = 3 (green 
dot-dashed line). Reionization photoevaporates gas in 
low-mass halos, steepening the PDF at log A = 1 and 
redistributing baryons back into the log A < 1 diffuse 
IGM reservoir. 


the volume and are able to keep a significant frac- 
tion of their hydrogen in neutral form. It is the 
LLS opacity that causes the mean free path of 
LyC radiation to remain relatively small even after 
overlap. At z < 6, the mean free path of ionizing 
photons in the IGM can be measured directly from 
the shape of the transmitted flux profile blueward 
of the Lyman limit in the stacked spectra of high- 
redshift quasars (Prochaska et al., 2009; Worseck 
et al., 2014; Becker et al., 2021). 


2.7 A clumpy IGM 


The IGM is highly inhomogeneous owing to the 
nonlinear evolution of density fluctuations over 
time. Inhomogeneities are key to understanding 
absorption features in the spectra of high-redshift 
sources, increase the effective recombination rate 
of the IGM, slow down the reionization process, 
and cause the mean free path of LyC photons 
to remain sub-horizon even after the epoch of 
overlap. A statistical description of the baryonic 
density field is provided by the volume-weighted 
probability density distribution (PDF) of IGM gas 
Py(A), normalized according to [5° Py(A)dA = 
1 (Miralda-Escudé et al., 2000). 

Figure 5 shows the PDF per unit log A com- 
puted from a large set of cosmological smoothed 


particle hydrodynamics (SPH) simulations in 
which a uniform UV background is switched on at 
z =9 (Bolton & Becker, 2009). The density PDF 
can be accurately parameterized by the linear rms 
density fluctuation and the local slope of the den- 
sity power spectrum at the scales of interest (Chen 
et al., 2022); that reference also provides fitting 
formulas that can be used in practical calculations. 

As shown by Pawlik et al. (2009), the presence 
of an ionizing background has a strong impact 
on the PDF. In the absence of photoheating, 
gravitational collapse proceeds unimpeded to the 
smallest scales set by the CMB decoupling, caus- 
ing a flattening of the PDF for overdensities 1 < 
log A S 2, a deficit of gas in the diffuse IGM 
(log A < 1), and shifting its maximum to lower 
overdensities. Reionization photoevaporates gas in 
low-mass halos, reducing the clustering at ~ 100 
kpe scales and steepening the PDF at logA > 1 
and redistributing baryons back into the log A < 1 
diffuse reservoir. The details are sensitive to the 
timing of reionization, the photoheating energy 
input per baryon, and radiative transfer effects like 
the self-shielding of gas inside halos. 

More detailed information on the clumpiness 
of the IGM can be extracted by second-order 
statistics like the 1D flux power spectrum Pr(k), 
which measures correlations of fluctuations in the 
transmitted flux fraction F' (the observed quan- 
tity) along the line of sight to distant quasars 
(e.g. Croft et al., 1998; McDonald et al., 2000; 
Irsié et al., 2017; Chabanier et al., 2019; Boera 
et al., 2019). These are caused by variations in 
the Lyman-a scattering optical depth 7 = —InF, 
and therefore probe the gas density, temperature, 
and velocity field of neutral hydrogen. After mea- 
suring the distance between pixels in units of the 
local velocity scale v and defining the flux con- 
trast Op(v) = F(v)/(F)—1, where (F’) is the mean 
transmitted flux, one can decompose each absorp- 
tion spectrum into Fourier modes p(k). Their 
variance as a function of the Fourier wavenum- 
ber k = 27/v is the flux power spectrum over the 
velocity interval Av, 


P(k) = Av(6r(k)’), (26) 
which is commonly expressed in terms of the 
dimensionless quantity AZ.(k) = kP(k)/7. Figure 


6 shows measurements of P(k) over the redshift 
range 4 < z < 5 from observations by eBOSS 
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(Chabanier et al., 2019), Keck, and the VIT tele- 
scopes (Irsi¢ et al., 2017; Boera et al., 2019). The 
figure also shows the best-fit P(k) from the cos- 
mic photoionization and photoheating history of 
Villasenor et al. (2022) (see also Fig. 3). The 
flux power spectrum contains information that are 
encoded across different velocity scales. While on 
scales larger than ~ 100 km s~! P(k) is sensi- 
tive to the ionization fraction of hydrogen in the 
IGM, on smaller scales the power spectrum is sup- 
pressed by the Doppler broadening of absorption 
features, the pressure smoothing of gas density 
fluctuations, the smearing effect of peculiar veloci- 
ties or, alternatively, by the free streaming of dark 
matter particles in the early Universe (e.g. Viel 
et al., 2013; Garzilli et al., 2019). 


3 Numerical Methods 


Numerical algorithms for modeling the dynam- 
ics of both collisionless dark matter particles and 
baryons are well established and mature, and we 
do not discuss them here (see Teyssier, 2015; 
Springel, 2016, for recent reviews that cover this 
topic in depth). We also do not survey models for 
star formation and stellar feedback, as these vary 
widely between different reionization simulation 
projects. Reionization simulations largely adopt 
recipes developed for and used in mainstream 
galaxy formation simulations. Recent reviews by 
Teyssier & Commergon (2019) and Vogelsberger 
et al. (2020) provide a general description of most 
of these recipes. 


3.1 Radiative transfer 


Radiative transfer algorithms solve for the evo- 
lution of the radiation field, taking into account 
emission, absorption, and scattering processes. In 
general, the evolution is described by the differ- 
ential equation (15) for the specific intensity I,, 
which is a function of seven variables: 3D posi- 
tion, 2D angular coordinates, time, and frequency. 
Because of the high dimensionality of the problem, 
radiative transfer is usually the most computa- 
tionally expensive component of a reionization 
simulation. In order to implement a radiative 
transfer scheme into cosmological simulations with 
a mass and spatial resolution comparable to those 
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Fig. 6 Dimensionless transmitted flux power spectra of the Lyman-a forest at redshifts z = 4— 5. The data points with 
error bars are from Irsié et al. (2017); Chabanier et al. (2019); Boera et al. (2019). Also shown are the maximum likelihood 
fit (black line) and 95% confidence interval (gray band) of Villasenor et al. (2022), computed from hundreds of cosmological 
hydrodynamics simulations with varying photoionization and thermal histories of the IGM. Note that the velocity scale 
corresponding to a given comoving spatial scale changes with redshift. The fractional differences between the observations 
and the best-fit model are shown in the bottom part of each panel. (Used by permission from B. Villasenor.) 


featured in modern galaxy formation runs, a radia- 
tive transfer algorithm should scale close to lin- 
early with the number of resolution elements, just 
like N-body and hydrodynamics algorithms. To 
satisfy this requirement, some degree of physical 
approximation and computational optimization is 
needed. 

Existing algorithms for solving the radiative 
transfer equation (15) can be broadly divided 
into two categories: “ray tracing” (we also include 
Monte-Carlo methods in this category) and 
“moments methods”. It is interesting to note that 
only the moments methods are being used in 
the largest fully-coupled simulations, and so it 
appears, at least superficially, that ray-tracing 
techniques may be gradually falling out of use. 


3.1.1 Ray tracing 


Ray tracing is a domain-specific jargon for the 
standard method of characteristics (e.g. Courant 
& Hilbert, 2008) for solving partial differential 
equations. The key advantage of using the method 
of characteristics in radiative transfer calculations 
is that photons move along geodesics (i.e. straight 
lines in the Newtonian limit); therefore one does 
not need to solve for the shape of characteristic 
curves, but only for the time dependence of the 
radiation intensity along each characteristic. 
Even with this simplification, ray-tracing tech- 
niques suffer from unfavorable scaling. In the 
simplest case of N, sources on a uniform grid of 
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size N?, and requiring that each grid cell has at 
least one ray from each source passing through 
it, the number of rays required is N, x 6 x N@. 
Given that each ray, on average, passes through 
N,/2 cells, the total number of operations then 
scales as N, x N?. In a large-box simulation, the 
number of sources may be substantial. Moreover, 
at fixed mass resolution, galaxies of a specified 
mass contain the same number of resolution ele- 
ments, and the number of sources N, at a given 
time is then a fixed (small) fraction of the total 
number of resolution elements. Therefore, for a set 
of simulations having fixed mass resolution but in 
progressively larger volumes, the cost of the sim- 
plest implementation of ray tracing scales as N° 
(the same scaling as a direct summation gravity 
solver), which is prohibitively expensive for any 
realistic simulation. 

A significant effort has gone into reducing 
this unfavorable scaling. Abel & Wandelt (2002) 
noticed that if 6N? rays are cast to cover all 
cells at the faces of the computational volume, 
multiple rays will pass through each cell close 
to the sources. They used the HEALPix subdivi- 
sion of the sphere (Gorski, 1999) to implement a 
more efficient adaptive ray-tracing scheme, where 
the angular resolution continually increases far- 
ther away from sources. For a given source, a small 
number of parent rays are cast and travel a short 
distance before they split into four daughter rays. 
Successive generations of splitting continue down- 
stream. This algorithm, however, still scales as 


N, x N3 with a smaller coefficient in front. Other 
variations of ray-tracing methods include group- 
ing near sources in a tree-like fashion (Razoumov 
& Cardall, 2005), merging closely aligned rays 
(Trac & Cen, 2007), or limiting the splitting of 
rays (McQuinn et al., 2007) in the Abel & Wan- 
delt (2002) scheme. All these approaches break 
the unfavorable N, x N? scaling at the expense 
of angular resolution. The effect of such resolution 
loss on the accuracy of the solution has not been 
thoroughly studied. 

There exist several variations of the ray-tracing 
method, e.g. “short characteristics” and “long 
characteristics”, which we do not cover here in 
detail as the differences between them are not 
fundamental (i.e. they do not eliminate the unfa- 
vorable scaling with N, x N?). The primary 
difference between the “short” and “long char- 
acteristics” approaches is that in the latter the 
rays originate from the source, while in the former 
the rays are centered on each target cell — such 
constrast is analogous to the difference between 
separate methods for solving hydrodynamics on 
the grid, like a Riemann solver versus a “Total 
Variation Diminishing” (TVD) scheme. 

A somewhat more distinct flavor of ray trac- 
ing is a class of Monte-Carlo methods. In this 
approach a source emits a number of “photon 
packets” in random directions, therefore randomly 
sampling a subset of all possible rays. As the 
simulation progresses, all possible directions get 
eventually sampled, but only at a finite tempo- 
ral cadence. These schemes, therefore, are less 
suitable for modeling rapidly varying sources. 


3.1.2 Moments methods 


The moments method converts the radiative 
transfer equation (15) into a hierarchy of moments 
of the radiation intensity I, 
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where J, is the mean specific intensity defined in 
equation (2), Fi = f dQn'I, is the photon flux, 
Ss is the angle-averaged radiation source function 
S, from Equation (15), and Q? is the flux source 
term. The second moment P?) = [dQn'n’I, is 
the radiation pressure tensor, and its trace is equal 
to Jy. The radiation Eddington tensor, defined as 


ij 
P, 


nid 
Vv Jy ? 


(28) 


is dimensionless and has a unit trace. In the gen- 
eral case the Eddington tensor or the radiation 
pressure tensor cannot be specified from the pho- 
ton energy density and flux alone. One can write 
down an equation for P’) similar to the above par- 
tial differential equations for J, and F%, but that 
equation would depend on the third moment of [,. 
The hierarchy of moments is therefore not closed. 

In direct analogy with gas dynamics, for scales 
much larger than the photon mean free path 
and when photons preferentially scatter rather 
than being absorbed, one can adopt a “diffu- 
sion” approximation under which h’? = 6° /3. In 
modeling reionization this limit is never achieved, 
however, and this approximation is not useful. 
Another interesting limiting case is when the gas 
is optically thin (k, = 0). Equations (27) then 
have an exact quadrature solution, which in the 
simplified case of isotropic sources (S, = 5, and 
Q!, = 0) takes the familiar 1/(47r?) inverse-square 
law form for the flux, 


Since the Eddington tensor cannot be computed 
exactly within the scope of the moments method, 
approximate schemes have been commonly used. 
Of these the two most widely used are the “M1 
closure” and “Optically Thin Variable Eddington 
Tensor” (OTVET) approximations. 

The M1 closure was first introduced by Lev- 
ermore (1984) with the following ansatz for the 


Eddington tensor (dropping the frequency depen- 
dence for brevity): 


l-a,.,, 38a-1 ; ,; 
aj 4 A) 
2 6 T 5) mn, 


hid = 


where n' = F'/||F|| is the unit vector in the 
direction of the flux propagation, 


3+4f? 


= pa faa 


and f = ||F||/J. In the vicinity of a single strong 
source f is close to 1, a@ is also close to 1, and 
hY = n'nJ, which is the correct solution in this 
case. In the opposite limit when multiple sources 
contribute more-or-less equally to the radiation 
field at a given location, then f < 1, a = 1/3, and 
h') x 5 /3 as in the diffusion limit. The functional 
form of a(f) then smoothly interpolates between 
the two regimes. 

In the OTVET approximation, first introduced 
by Gnedin & Abel (2001), the Eddington tensor 
is computed as the ratio of the optically thin radi- 
ation pressure PJ and mean intensity J, from 
equations (29). OTVET then uses this “optically 
thin” tensor in the full, “optically thick” equations 
for the radiation mean intensity and flux. In the 
vicinity of a strong source (or anywhere in the 
simulated volume in the case of a single source), 
OTVET produces the same (exact) answer as the 
M1 closure, h’) = n‘n’. In a general location with 
several sources of radiation, OTVET accounts for 
their spatial location and clustering, but it may do 
so incorrectly in the presence of luminous sources 
embedded in very optically thick regions — in that 
case OTVET would still include their contribution 
to the Eddington tensor. 

It is important to keep in mind that Eqs. (27) 
are written in conservative form: in the absence 
of absorption the radiation density and flux are 
conserved, and in the presence of absorptions they 
are reduced proportionately. Thus, errors in the 
Eddington tensor do not violate photon density or 
flux conservation, but rather result in the photon 
flux being advected in the wrong direction. For 
example, if the recombination term in equation 
(33) is not large, the overall reionization history 
is not affected by errors in the Eddington ten- 
sor (provided that the LyC mean free path is 
short, as it is the case during reionization) since 
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the number of ionizations is captured correctly 
by equations (27). For the same reason these 
equations automatically yield the correct speeds 
for expanding ionization fronts. In the presence 
of large gas density fluctuations and significant 
radiative recombinations, the different approxima- 
tions for computing the Eddington tensor may 
produce somewhat different results, as the pho- 
ton flux may be incorrectly advected into dense 
regions in one scheme (where it is absorbed by 
gas that later recombines) or into lower density 
regions in the other, where it results in per- 
sistent photoionizations. A detailed comparison 
between M1 closure and OTVET has never been 
performed, however. 

An interesting variation of the moments meth- 
ods was developed by Finlator et al. (2009), who 
used ray tracing to compute the Eddington ten- 
sor accurately, without any approximation. This 
approach, however, combines the computational 
expense of both the ray-tracing and the moments 
schemes. For comparison, Finlator et al. (2018) 
carried out radiative transfer simulations on a 
uniform 64° grid, while modern reionization simu- 
lations with the moments methods, such as CROC 
and THESAN, typically use 2000° resolution ele- 
ments and achieve a dynamic range with AMR or 
moving mesh in excess of 10°. 


3.1.3 Reduced speed of light 
approximation 


So far we have only discussed techniques for 
addressing the spatial dependence of the radiative 
transfer equation. Its time dependence presents 
an additional challenge, since photons travel at 
the speed of light, which is much larger than any 
other relevant speed in a typical reionization sim- 
ulation. In any explicit scheme, the time-step will 
be limited by an appropriate Courant—Friedrichs— 
Lewy stability condition (Courant et al., 1928) 
to a value that is about Umax/c times smaller 
than the hydrodynamic time-step (here Umax is the 
maximum gas velocity in the simulation, includ- 
ing the sound speed). For reionization simulations 
this ratio can be as small as 10~%. Clearly this 
is challenging, as it is computationally expensive 
(although not impossible) to make 1,000 radiative 
transfer time steps for every hydrodynamic time 
step in a realistic simulation. 


The time dependent term in the radiative 
transfer equation (15) is only crucial for capturing 
“light fronts” — waves in the radiation field that 
propagate with the speed of light just like waves in 
gas propagate with the speed of sound. Presently 
it is unclear if capturing such waves is ever impor- 
tant.” For example, ionization fronts in the general 
IGM typically propagate at speeds of order a few 
thousand km/s. Capturing these ionization fronts 
is sufficient for following the overall distribution 
of ionized bubbles, and if that is the primary goal 
of a reionization simulation, capturing light fronts 
may be an overkill. Following this logic, Gnedin 
& Abel (2001) introduced the “Reduced Speed of 
Light” (RSL) approximation, which replaces the 
speed of light c in the first term of equation (15) 
with an “effective” speed of light ¢, 


1d0J, x lod, 
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For example, taking ¢ = c/10 would still be suf- 
ficient for capturing ionization fronts propagating 
at 3,000km/s with 10% precision, while increas- 
ing the time-step in the radiative transfer solver 
by a factor of 10. Obviously, such a modification 
would not be appropriate when light fronts need 
to be tracked or when I-fronts move at a signifi- 
cant fraction of the speed of light. However, this 
is not expected to be the case during reionization 
(Shapiro et al., 2006). 

Unfortunately, the RSL approximation is a 
delicate numerical tool, and there exist multiple 
examples of its incorrect use. Reducing the speed 
of light in the transfer equation for the specific 
intensity of the cosmic UV background results in 
the incorrect redshifting of photon frequencies, so 
the RSL approximation can only be used when 
computing the contribution from sources located 
at distances smaller than the cosmic horizon. The 
radiation field must then be represented as the 
sum of two components, one from local sources 
and the other from “distant” ones — sources 
outside the simulation box. In simulations that 
model the radiation field with just one spatial 
component, the RSL approximation should never 
be used! The second pitfall is that only the speed 
of light appearing in the first term of equation (15) 


?Note that radiation sources can also change rapidly. This, 
however, is presumably captured by the simulation with, for 
example, a hydrodynamic solver or a star formation algorithm. 
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can be reduced. The speed of light enters explic- 
itly or implicitly in many physical expressions, and 
thus great care must be taken in deciding where it 
can and where it cannot be modified. A detailed 
description of the RSL approximation and how it 
should be implemented can be found in Gnedin 
(2016b). 


3.2 Thermodynamics 


The cooling and heating processes that control the 
thermodynamics of cosmic gas are well known, but 
complex. For simple geometries they can be mod- 
eled in great detail with the “industry-standard” 
Cloudy software (Ferland et al., 1998). In realistic 
simulations, however, it is still too computation- 
ally expensive to run a code like Cloudy for every 
resolution element at every step, and hence one 
has to resort to some reasonable approximations. 

The complexity of the problem is highlighted 
by considering the dependence of the cooling 
function A and the heating function H in the 
temperature equation (6). In addition to obvious 
arguments such as gas temperature T and total 
hydrogen density ny, in the most general case 
the cooling and heating functions depend on a 
large number of other parameters that describe 
the fractional abundance for each atomic, ionic, 
and molecular species at each of their quantum 
level. If all species can be assumed to be in local 
thermodynamic equilibrium (almost always a safe 
assumption) and to be close enough to ionization 
equilibrium (sometimes a safe assumptions, see 
discussion below), and the abundance pattern for 
heavy elements is assumed to be the same every- 
where in the simulation volume (e.g., to be solar), 
then the cooling and heating functions become 
just functions of T’', ny, the gas metallicity Z, and 
the radiation field J,, 


{H, A} = {H,A}(T, nu, Z, Jv). (30) 


Even under these simplifying assumptions 
equation (30) is too complex, since its RHS is not 
a function but rather a functional, as it depends 
on the specific mean intensity J). 

One approach to account for this dependence 
is to use a spatially uniform radiation spectrum. 
For example, the spectrum of the evolving cos- 
mic radiation background has been used in several 
galaxy formation simulations (Kravtsov et al., 


2002; Wiersma et al., 2009; Vogelsberger et al., 
2013; Vogelsberger et al., 2014; Schaye et al., 
2014) — so the dependence of the heating and 
cooling functions on the radiation spectrum is 
reduced to a dependence on cosmic time only. The 
main limitation of such an approach is that the 
resultant cooling and heating functions are only 
valid in the general IGM, but become inaccurate 
near and inside galaxies, where accounting for gas 
thermodynamics is especially important. 

An alternative simplification is to only treat 
hydrogen and helium accurately. In that limit 
not only the cooling and heating rates can be 
computed accurately, but even the assumption of 
ionization equilibrium can be relaxed — for mod- 
eling reionization this is particularly relevant, as 
equilibrium is not maintained at ionization fronts 
(Gnedin, 2000; Aubert & Teyssier, 2008; Gnedin, 
2014; Kannan et al., 2021). The metal-dependent 
part of the cooling and heating functions is then 
added to the cooling and heating processes from 
hydrogen and helium with some approximation 
— for example, using a uniform cosmic radiation 
background or an approximation that parameter- 
izes the effect of an arbitrary radiation field with 
several photoionization rates (Gnedin & Hollon, 
2012). 

The assumption of ionization equilibrium for 
all metal ions may not hold in all relevant phys- 
ical regimes. Richings & Schaye (2016) explored 
departures from ionization equilibrium in galaxy 
formation simulations with a full time-dependent 
modeling of multiple metallic species, and found 
some differences. It is currently unknown if these 
differences have any systematic effects on gas 
dynamics on large scales. 


4 A General Overview of the 
Reionization Process 


The reionization of intergalactic hydrogen during 
the first billion years of cosmic history involves 
the early production, propagation, and absorption 
of LyC radiation. It is, in essence, a relatively 
simple story of UV sources and sinks. While it 
is generally agreed that star-forming galaxies and 
AGNs are the best candidates for providing the 
bulk of the photoionizing power across cosmic 
time, the apparent contrast between the sharp 
decline in the number density of bright quasars at 
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z 2 3 (e.g. Kulkarni et al., 2019) and the relatively 
shallower decline in the UV luminosity density 
from galaxies (e.g. Madau & Dickinson, 2014) 
has led to the predominant theory that the bulk 
of the photon budget for hydrogen reionization 
came from massive stars. There also remains some 
room for a contribution from more exotic possi- 
bilities such as dark matter annihilation or decay 
(Mapelli et al., 2006; Liu et al., 2016), primordial 
globular clusters (Ricotti, 2002; Ma et al., 2021), 
mini- (Madau et al., 2004) and micro-quasars 
(Mirabel et al., 2011). The basic energy source 
in all these scenarios is often rather different — 
rest-mass and nuclear binding energy in the case 
of particle annihilation/decay and stellar nucle- 
osynthesis, and gravitational binding energy in the 
case of accretion onto black holes. And various 
reionization sources will be characterized by very 
distinct radiative efficiencies, spectra, abundances, 
and spatial distributions. 

The sinks of LyC photons during reionization 
come essentially in two flavors: hydrogen atoms 
in the diffuse IGM (both pristine and previously 
photoionized and recombined) and the LLSs (we 
use this term here to also include the rarer and 
denser HI Damped Lyman-a absorbers). The sep- 
aration between the two is, admittedly, somewhat 
artificial, as there is just one continuous HI den- 
sity field in the Universe. Nevertheless, the LLSs 
are well-defined observationally (they are caused 
by structures with integrated HI column densities 
large enough to produce optically thick absorp- 
tion at the Lyman edge) after the reionization of 
the diffuse IGM is completed, and therefore one 
can extrapolate their abundance into the reioniza- 
tion epoch and draw a distinction between these 
absorbers and hydrogen atoms closer to the mean 
density. The key physical difference is that LLSs 
trace high density regions in the outskirts of galax- 
ies (Fumagalli et al., 2011; Shen et al., 2013) that 
occupy a small portion of the volume and are able 
to keep a significant fraction of their hydrogen in 
neutral form for a Hubble time. By contrast, a 
comparable amount of hydrogen gas in the diffuse 
IGM is photoionized almost instantly.° 


3 For example, a region of physical size 100 kpc at the mean 
density contains as much gas as a 10 kpc LLS with overdensity 
A = 1,000. A typical ionization front moving with 3,000 km 


s would ionize such diffuse region in 30 Myr. 


Fig. 7 A cartoon illustration of the relationship between 
the sizes of HIT bubbles and the LLSs they contain. The 
green color marks neutral gas, while ionized gas is in black. 
The small ionized bubble on the left keeps expanding, as 
many LyC photons (yellow rays, emitted from the source 
at the center) are absorbed by the surrounding neutral gas 
rather than by the few LLSs it contains. The large bubble 
on the right has instead stopped expanding, as it contains 
many LSSs that determine the mean free path of LyC 
radiation. 


On quite general grounds, one expects that a 
fraction of the ionizing radiation emitted by indi- 
vidual sources will leak into intergalactic space 
and carve out cosmological HII regions (more 
commonly colloquially called “ionized bubbles” ) 
in their immediate vicinity. These isolated ionized 
bubbles will expand in an otherwise largely neu- 
tral phase as more LyC photons are produced, and 
will start overlapping — in a “swiss-cheese” reion- 
ization topology — when their typical size becomes 
comparable to the mean distance between highly 
clustered groups of sources; for galaxies this would 
be comparable to their correlation length, about 
5 Mpc (Barone-Nugent et al., 2014). Merged bub- 
bles will contain more radiation sources, and the 
ionization fronts delimiting their boundaries will 
expand faster as a consequence. A faster expand- 
ing bubble will rapidly engulf more sources and 
expand even faster, in a runaway process that con- 
tinues until the mean free path of LyC photons 
becomes effectively controlled by the LLSs rather 
than by the typical size of HII bubbles (see Fig. 7). 
There must therefore be a short time interval when 
the mean free path increases sharply (and so does 
the UV background intensity) causing a sudden 
drop in the neutral fraction. This unambiguous 
“epoch of overlap” can then be defined mathe- 
matically as the moment when the LyC mean free 
path increases at the fastest rate (Gnedin, 2004), a 
physically-motivated definition that is preferable 
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Fig. 8 The redshift-dependent Amfp of H-ionizing radia- 
tion measured from quasar spectra. Figure 8 from Fan et al. 
(2006) is overplotted with more recent observational deter- 
minations from Worseck et al. (2014, green) and Becker 
et al. (2021, blue). The plot from Fan et al. (2006) also 
shows the evolution of the mean free path in the early sim- 
ulations of Gnedin (2004). Note how the hightest-redshift 
value lie well below extrapolations from lower redshifts, a 
rapid decrease in the mean free path that is qualitatively 
consistent with models wherein reionization ends at z ~ 6. 
In the units on the y—axis, the symbol h denotes the Hub- 
ble constant in units of 100 km s~t Mpc~t. (Used by 
permission from Xiaohui Fan.) 


to some arbitrary chosen thresholds for the mass- 
or volume-weighted hydrogen neutral fractions. 
Under some assumptions, the LyC mean free 
path can be inferred from the absorption spectra 
of high redshift quasars. Two recent measure- 
ments are shown in Fig. 8 together with the first 
determinations by Fan et al. (2006). While more 
observations are needed to fully constrain the evo- 
lution of the IGM opacity at these redshifts, the 
existing data already hint at the rapid decrease in 
the mean free path expected for a epoch of over- 
lap at z ~ 6. There is one exception to the above 
scenario, and that is when the number density of 
LLSs is so high that the growth of HII regions 
stops before overlap and overlap never actually 
occurs. The process of subsequent reionization 
(and whether it completes at all) is then controlled 
entirely by the evolution of the LLS abundance. 
As the observational data in Fig. 8 indicate, this 
situation is not realized in our Universe. 


5 Analytical Models of 
Reionization 


To zeroth order, tracing the history of cosmic 
reionization is mainly a photon counting exercise, 
where the growth rate of HII regions is equal 
to the rate at which ionizing photons are pro- 
duced minus the rate at which they are consumed 
by radiative recombinations. Most of the existing 
analytical tools, while appearing rather dissimi- 
lar from each other, apply in essence the same 
bookkeeping algorithm to different spatial scales. 
For example, one of the long running historical 
debates in the topology of reionization is whether 
such process proceeds “inside-out” (high density 
regions are ionized first, and lower-density regions 
follow, Furlanetto et al. 2004) or “outside-in” 
(lower density regions are ionized first, Miralda- 
Escudé et al. 2000). While the two approaches 
have been often presented as mutually exclusive, 
we shall show in this section that they are, in fact, 
just two different stages of the same “inside-out- 
back in” process. 


5.1 The “reionization equation” 


Let us start then with counting the number of 
hydrogen ionizing photons present in an arbitrary 
comoving patch of the Universe (without regard 
to its size, shape, gas density, or other physical 
properties). This patch is fully ionized when the 
number of photons with energies above 13.6 eV, 
N,, either emitted inside the region or arriving 
into it from external sources, is equal or exceeds 
the number of hydrogen atoms Ny in the patch, 
Ny > Nu. If this condition does not hold the patch 
is only partially ionized, with Nyy < Ny and 


Ny = Nut. 


This is not enough, however, since atoms can 
recombine following the absorption of UV radi- 
ation, and must be photoionized again. Hence, 
our “master equation” should also include photon 
losses from radiative recombinations, 


dt! 
trec (t’) , 


N,(t) = Nyr(t) + [ Nuit’) (31) 
0 


where tye. is the recombination time inside our 
patch. Reionization is complete when Nyy = Ny 
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in every region of the Universe that is close (in 
some particular sense) to the mean cosmic density. 
Averaging now our master equation over a repre- 
sentative large comoving volume of the Universe 
V, denoting this average with angle brackets, and 
differentiating with respect to time, one gets 


ld 1d 1 / Nan 
—_— — eee . 2 
(N,) V ag Num) + = ( E ) (32) 


For historical reasons, the universe-average ion- 
ized fraction is often labeled with the symbol 


Q, 


N 
Q = Tn HII = (nu I) ; 
Vo Ny (nu) 


Rearranging terms, one recovers the original 
“reionization equation” of Madau et al. (1999): 


d ni 
Q = 10n Qo (33) 
dt TH,0 trec 

where tree is an “effective”? recombination 


timescale and ny,o = 1.9 1077 cm~? is the mean 
comoving cosmic density of hydrogen atoms. Here, 
we have written the comoving volume-averaged 
emission rate into the IGM of ionizing photons as 


<p hes 
Nion = yy) = fesc&ionPSFR- (34) 


In reionization studies, it is convention to include 
in the ionization balance equation only a frac- 
tion of the hydrogen photoionizations and radia- 
tive recombinations, i.e. those that take place 
in the diffuse, low-density IGM. Photons that 
are absorbed locally, within the immediate high- 
density vicinity (the ISM and the CGM) of ion- 
izing sources, are accounted for by a reduction of 
the source term via the escape fraction parame- 
ter fese (see also Section 2.5). The third equality 
in Equation (34) assumes a galaxy-dominated ion- 
izing emissivity, with pspr(t) being the cosmic 
star formation rate density, €;., the appropriately- 
averaged number of LyC photons emitted per unit 
SFR (a “yield” determined by stellar physics and 
the shape of the initial mass function), and fog. the 
fraction of such photons that leaks into the IGM. 

Some care must also be exercised with the 
photon sink term in Equation (33). Since the 
recombination rate is quadratic in density, it 
depends on the actual gas density distribution 


within the volume and cannot be expressed via 
globally averaged quantities only. It is custom- 
ary then to introduce a “recombination clumping 
factor” Cy and write 


1 n nad 
f= (MHF) _ (aft) (ne)Cun, (35) 
trac (nu II) 
where ( A ) 
_ NA UINe Ay TY 
Cun = ae 
(nu) (ne) (aft) 


is in general time-dependent. Neglecting density- 
dependent temperature fluctuations, which can 
impact the recombination coefficient, and assum- 
ing Ne K NYT, ONe can rewrite the clumping factor 
in the more familiar form Cy yy = (n%71)/(nu)?- 

The clumping factor is an external parameter 
that must be computed with the help of cosmo- 
logical hydrodynamics simulations (Kohler et al., 
2007; Pawlik et al., 2009; Finlator et al., 2012; 
Shull et al., 2012b; Kaurov & Gnedin, 2014; So 
et al., 2014). It has to be chosen appropriately to 
avoid double counting, since the impact of recom- 
binations within the host halos of ionizing sources 
is already implicitly incorporated as a reduction in 
the source term through the escape fraction fose. 
For example, all volume elements denser than a 
given gas overdensity threshold — say A > 100 — 
should be removed from calculations of the clump- 
ing factor if ionizing photons were counted as 
escaped into the IGM once they entered regions 
of gas overdensities A < 100. This means that 
the escape fraction and the clumping factor are 
not independent, but as long as the definitions 
of escape fraction and clumping factor refer to 
the same decomposition of the gas into IGM and 
ISM/CGM, the value of the overdensity threshold 
separating these different components can be cho- 
sen arbitrarily. Note, however, that the exclusion 
in the ionization balance equation of gas above a 
certain density implies that the quantity nyo in 
Equation (33) should not be, strictly speaking, the 
mean hydrogen density of the Universe but rather 
the redshift-dependent mean density of its diffuse 
component. 

Equation (33) above statistically describes the 
transition from a neutral Universe to a fully 
ionized one. Extensively used in the literature 
(e.g., Haardt & Madau, 2012; Kuhlen & Faucher- 
Giguere, 2012; Bouwens et al., 2015; Robertson 
et al., 2015; Khaire et al., 2016; Ishigaki et al., 
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2018; Sharma et al., 2017) as it allows an esti- 
mate of the photon budget required to achieve 
reionization with a fast exploration of parameter 
space, this simple ODE has been shown to provide 
an acceptable description of the results of radia- 
tive transfer simulations (Gnedin, 2000, 2016a). 
The main limitation of this approach is that it 
assumes that all LyC photons escaping from indi- 
vidual galaxies are absorbed by the diffuse IGM, 
mathematically permitting unphysical values of 
Q above unity when reionization is completed 
and the Universe becomes transparent, because 
our master equation (31) is only valid until the 
patch is fully ionized.* Another limitation of the 
reionization equation is that it ignores collisional 
ionizations. As long as these take place inside 
virialized halos, they can be accounted for by 
appropriately choosing the parameter fesc; how- 
ever, collisional ionizations in filaments outside 
virialized halos are not included in this approach. 

To provide an illustrative description of how 
the global reionization history of the Universe may 
proceed over cosmic time, we have numerically 
integrated Equation (33) from z = 10 onwards, 
assuming a constant emission rate of LyC photons 
per hydrogen atom of nion/nu,o = 2.7 Gyr. In 
this toy model, the IGM is defined as all gas at 
an overdensity threshold of A < 100. Its effective 
recombination timescale is 


tae = 2.3((1 + 2)/6]* Cyt, (37) 


a fitting formula provided by So et al. (2014) 
and based on the analysis of a radiation hydro- 
dynamical simulation of reionization that begins 
at z = 10 and completes at z ~ 5.8. The result- 
ing average ionized fraction, shown in Figure 9, 
appears to be consistent with a number of obser- 
vational constraints on the state of the z > 5IGM, 
and shows that the process of reionization is quite 
extended: even if it completes at redshift < 6, 
40% of hydrogen atoms in the Universe are ionized 
already at z ~ 7.5. Note that we have artificially 
set Q = 1 when reionization is completed, the 
neutral fraction in the IGM drops abruptly, and 


“See Madau (2017), however, for a revised formulation where 
the source term in Equation (33) is modified to explicitly 
account for the presence of the optically thick LLSs that trace 
the CGM (Crighton et al., 2019) and determine the mean free 
path of ionizing radiation after overlap. 
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Fig. 9 A toy model of reionization history obtained by 
integrating Equation (33) with a constant emission rate 
of ionizing photons per hydrogen atom, nion/ny,9 = 2.7 
Gyr—!. The average ionized fraction Q is shown for an 
IGM that recombines according to Equation (37)). The 
data points represent constraints from the dark pixel statis- 
tics at z = 5.6 and z = 5.9 (blue squares, McGreer et al. 
2015), the gap/peak statistics at z = 6.32 (magenta square, 
Gallerani et al. 2008), the damping wing absorption pro- 
files in the spectra of quasars at z = 6.3 (orange pentagon, 
Schroeder et al. 2013) z = 7.1 (turquoise hexagon, Greig 
et al. 2017; Mortlock et al. 2011), z = 7.29 (red triangle, 
Greig et al. 2022), and z = 7.54 (green circle, Bafiados et al. 
2018), the redshift-dependent prevalence of Lyman-a emit- 
ters (LAEs) at 7 < z < 8 (firebrick pentagons, Schenker 
et al. 2014, and purple square, Mason et al. 2019a), and 
their clustering properties at z = 6.6 (gold hexagon, Ouchi 
et al. 2010). 


photoionizations come into balance with radiative 
recombinations. 


5.2 PDF-based models 
(“outside-in” ) 


A different treatment was presented by Miralda- 
Escudé et al. (2000), who used the volume- 
weighted PDF of IGM gas Py(A) already dis- 
cussed in Section 2.7 to describe the ioniza- 
tion state of an inhomogeneous Universe. They 
assumed that inside each ionized bubble (and for 
the Universe as a whole after all the low density 
IGM is ionized) reionization proceeds gradually 
and “outside-in”, first into the underdense voids 
and then more gradually into overdense regions. 
At any given epoch, all the gas with density 
above a critical overdensity A,(t) remains neu- 
tral, self-shielded from the ionizing background, 
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while all lower density material is completely ion- 
ized. This ansatz cannot clearly be exact, since 
the degree of ionization of gas at a given density 
will depend on the environment, but it may pro- 
vide a useful approximation for understanding the 
reionization history of cosmic structure. The over- 
lap of HII regions occurs then first through the 
lowest density “tunnels” found between sources, 
both because fewer atoms need to be ionized per 
unit volume and because fewer recombinations 
take place, while the denser, neutral gas with A > 
A,(t) limits the mean free path of LyC radiation 
to sub-horizon scales. 

Given a PDF, the mass-averaged ionized frac- 
tion is 


Ai 
me | APy(A)dA. (38) 


The comoving volume-averaged recombination 
rate can then be written as 


A; 
aod nnenm) = 22 | A?Py(A)dA, (39) 
rec,0 JO 


since only ionized low-density regions contribute 
to the sink term. In the last expression, trec,o is 
the recombination time at the cosmic mean den- 
sity and we have again neglected the effect of 
the temperature-density relation on the recombi- 
nation coefficient. In this formalism, the clumping 
factor is simply 


Ai 
Cun = Fy’ | A?Py(A)dA. (40) 
0 


The critical overdensity A; for a particular ion- 
ized patch of volume V (which may be different 
in different patches) is then determined implicitly 
by our master equation (31), 


N,(t) = Fu (Ai(t))+ i 2 / A2Py(A,t)dA. 
) 0 

(41) 
Equation (3) from Miralda-Escudé et al. (2000) 
is then obtained by differentiating this equation 
with respect to time and dividing it by V and 
by the Hubble parameter H. For typical Py (A), 
the second term in Equation (41) dominates over 
the first term for sufficiently large A; (typically 


A; 2 10 at z ~ 6). Hence, this class of models 
is commonly called “outside-in”, with “out” and 
“in” referring to low and high density gas, respec- 
tively. The former gets ionized first, while the 
latter (excluding the ISM in the very vicinity of 
UV sources) is ionized later, as N, and A, increase 
with time. In the late stages of reionization, when 
N, is sufficiently large, it is recombinations in 
high density regions that consume most of the 
LyC photons and determine the rate at which 
reionization proceeds. 

There are two main limitations of the PDF- 
based approach: 1) the gas density distribution is 
a fit to the results of cosmological hydrodynamical 
simulations that are sensitive to resolution, box 
size, as well as thermal and photoionization his- 
tory. The PDF choosen by Miralda-Escudé et al. 
(2000), for example, 


(A-?/3 — Go)? 


Py (A) = Aexp | 20/3) 


| A-®, (42) 


was derived assuming that the initial density fluc- 
tuations form a Gaussian random field, the gas 
in voids is expanding at a constant velocity, and 
the density field is smoothed on the Jeans scale 
of the photoionized gas. Here, A,Cpo,d9 and £ 
are redshift-dependent best-fit parameters. This 
PDF has been shown to be inadequate at describ- 
ing the results of more modern simulations for 
log A > 1 (Pawlik et al., 2009; Bolton & Becker, 
2009); and 2) only recombinations leading to the 
consumption of LyC photons that escaped from 
the high-density ISM of star-forming regions con- 
tribute to the balance in Equation (41). The 
number of ionizing photons leaking from the ISM 
must still be parameterized by an escape fraction 
fese Which, contrary to the clumping factor, is not 
a derived quantity. Hence, the “outside-in” model 
is incomplete, and should be supplemented by an 
equation for fesc(Ai). 

PDF-based models have found wide applica- 
bility in studies of the post-reionization IGM as 
probed by observations of hydrogen absorption 
in the spectra of luminous high-redshift quasars. 
Ignoring the effects of peculiar velocities and of 
thermal broadening, the Gunn-Peterson (Gunn & 
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Peterson, 1965) optical depth to Lyman-a scatter- 
ing at redshift z is given by 


2 
Te 
T= mecH(z)/20anal (43) 


where f, is the oscillator strength of the Lyman- 
qa transition, \. = 1216A, H(z) is the Hubble 
parameter, and nyy is the proper local density of 
neutral hydrogen in the IGM. In photoionization 
equilibrium, the neutral hydrogen fraction and 
therefore the optical depth depend on the local 
density of the IGM. For a highly ionized region 
with overdensity A in photoionization equilibrium 
one derives 


(A) ~9x 10-4 (14 yx) 


where ['_j9 is the hydrogen photoionization rate 
in units of 10-!? s-!, F(z) = [Qyv(1 + 2)? + 
Qq]'/2, x accounts for the extra electrons liber- 
ated by helium reionization, and the temperature 
is approximately related to the gas density via 
a power-law density-temperature relation of the 
form T = TyA?—!, which should be valid for 
low-density gas (Hui & Gnedin, 1997). The mean 
normalized transmitted flux through the IGM 
(which is the observed quantity) is then 


(F) =(e7) = iy e TA) DP, (A)dA. (45) 


Thus, as long as the PDF and the temperature- 
density relation are sufficiently well understood 
from theory and simulations, the observed mean 
transmitted flux can be used to determine the 
hydrogen photoionization rate Py; (McDonald & 
Miralda-Escudé, 2001; Fan et al., 2002). Note that, 
at high redshift, the transmitted flux distribu- 
tion e~74) Py(A)A peaks in underdense regions 
with A ~ 0.3, as denser gas produces saturated 
Lyman-a absorption with zero transmitted flux. 
For illustrative purposes, we have calculated 
here the mean transmitted flux as a function of 
redshift using the polynomial fit to the PDF from 
Bolton & Becker (2009), and assuming constant 
Ty = 10° K,.9 =. 1.8, and Ty = 10-1 s—" (et 
Fan et al., 2006; Becker & Bolton, 2013; Gaikwad 
et al., 2017; Boera et al., 2019; Walther et al., 
2019; Gaikwad et al., 2020; Villasenor et al., 2022). 
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Fig. 10 Redshift evolution of the observed mean transmit- 
ted flux from Becker & Bolton (2013) (blue dots), Gaikwad 
et al. (2021) (turquoise pentagons), Bosman et al. (2018) 
(magenta squares) and (Eilers & Hennawi, 2018) (yellow 
pentagons). The dashed curve shows the prediction of a 
PDF-based model with constant Ty = 104K, y = 1.3, 
and Ty = 10—-12s~1 (see text for details). The data sug- 
gest a rapid change in the ionization state of the IGM at 
z > 5 produced by a decrease in the ionizing ultraviolet 
background, an indication of the end stages of hydrogen 
reionization. 


Figure 10 compares the effective Lyman-a opti- 
cal depth, —In(F), predicted in the redshift range 
2 S$ z S 6 with the observations. Clearly, a con- 
stant hydrogen photoionization rate is unable to 
reproduce the rapid increase of the IGM opacity 
observed at z > 5, suggesting a swift transition in 
the ionization state of the Universe. Both the evo- 
lution of the observed Lyman-a opacity as well as 
the increased scatter of the measurements (Becker 
et al., 2015) are the most compelling constraints 
to date on the end of the epoch of reionization. 


5.3 Excursion set-based models 
(“inside-out” ) 


While able to characterize an inhomogeneous 
reionization process where overdense regions are 
photoionized after underdense voids, PDF-based 
models implicitly assume a uniform UV _ back- 
ground throughout the Universe and do not take 
into account large-scale fluctuations in the density 
field. The success of the extended Press-Schechter 
formalism for describing the correlation of peaks 
in the primordial mass distribution (Press & 
Schechter, 1974; Bond et al., 1991; Lacey & Cole, 
1993) has prompted efforts to develop a similar 
approach for studies of the topology of reioniza- 
tion on large scales. Pioneered by Furlanetto et al. 
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(2004) and Furlanetto & Oh (2005), such modeling 
has been subsequently used in a number of appli- 
cations, most notably in the 2lcmFAST code? for 
estimating the 21cm signal from neutral hydro- 
gen during the epoch of reionization (Mesinger & 
Furlanetto, 2007; Mesinger et al., 2011). The key 
idea here is to provide a statistical representation 
of the fluctuating ionization pattern on the scales 
that are most relevant to observations by connect- 
ing the number of ionizing photons present in a 
spherical region of comoving radius R or, equiv- 
alently, linear mass scale M = (47/3)pR? (where 
p is the mean density of the Universe) to the col- 
lapsed mass fraction f.oy in that region. According 
to the extended Press-Schechter formalism, this is 
given by 


Oc = dr(z) 


[Osim —oOR (z)?] 


eon = erfc (46) 


for Gaussian fluctuations on the scale R, where 
dp(z) and opr(z) are the linear overdensity in a 
region of size R and the rms density fluctuation 
in all regions of size R respectively, 6. is the 
linear-theory density threshold for collapse (equal 
to ~ 1.686 for spherical collapse), and Oyin is 
the rms density fluctuation on the smallest scale 
that can harbor an ionized source. It can be inter- 
preted as the spatial scale Rin corresponding to 
the mass scale My,i, for the minimum mass of an 
ionizing source, which, in turn, may correspond to 
atomic-cooling halos of virial temperature 10* K.° 
Furlanetto et al. (2004) postulated that the num- 
ber of hydrogen ionizing photons per hydrogen 
atom in the region is 


N. 


Na = Cfeolls (47) 


where the factor ¢ depends on uncertain source 
properties like the efficiency of LyC photon pro- 
duction, the escape fraction of these photons from 
the host galaxy, and the star formation efficiency, 
and can be a function of mass M. Then, from our 
master equation (31), the region under investiga- 
tion is fully “self-ionized” when it contains enough 


Shttps://github.com/21cmFAST/21cmFAST. 

SA pre-factor of 2 commonly added to this equation to 
achieve fio) 4 1 at t + oo can always be incorporated into 
the ¢ factor of Equation (47). 
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Fig. 11 Left panel: A thin slice through a 120 comoving Mpc simulation from the CROC suite. Each pixel in the image is 
color-coded by the value of redshift at which gas at that spatial location was 90% ionized, with the color-bar giving the actual 
values for each color. The inset zooms on a region that contains both early ionized bubbles and high density gas ionized 
later. Right top: Trajectories for the simulation box on the left at z = 8; blue lines (one solid, the other highly translucent) 
mark trajectories that are not yet ionized, read lines indicate ionized trajectories. The green line shows the reionization 
barrier. Right bottom: Distribution of initial overdensities on scales R/Lpox = 0.03. Again blue indicate neutral gas and 
red ionized; purple is the sum of the two. Both distributions are broad, so the barrier is “fuzzy”. One plausible definition 
of a reionization barrier is the overdensity at which the PDFs of neutral and ionized trajectories cross each other (vertical 


green line). 


mass in luminous sources, i.e. when 


dt 
Cfeou = 1 +f LHI: (48) 
0 


rec 


This, together with equation (46), places a con- 
dition dr > 6,(R,z) on the mean overdensity 
within an HII region of size R, where 06,(R, z) 
is the threshold or “ionization barrier”. At any 
given (randomly chosen) position in space the 
evolution of dz as a function of R resembles a ran- 
dom trajectory. Furlanetto et al. (2004) used the 
first-crossing distribution of random walks with a 
“moving” (because it is a function of the spatial 
scale R) barrier to model the size distribution of 
HII bubbles during reionization. 

A practical implementation of the barrier- 
crossing technique is as follows. Consider a spa- 
tially variable field of density fluctuations 6(x) (it 
can be either a linear Gaussian field or an actual 
nonlinear density field) at a given time in cosmic 
history, and an arbitrary position xg. Let dp(x) be 
the density field smoothed with a particular spa- 
tial filter of size R (the specific shape of the spatial 
filter is usually not very important). A sequence of 
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values dp(Xo) as a function of R forms a “trajec- 
tory” in the dg — R plane. The largest value of R 
for which 6R(xo) = 5p(R) is called the “first bar- 
rier crossing”. The Furlanetto et al. (2004) model 
then considers the sphere of size R centered on 
xo as fully ionized, since the number of ionizing 
photons per baryon inside that sphere exceeds the 
threshold value in Equation (47). A location for 
which dr < 6,(R) for any R is considered to be 
neutral at the cosmic time under consideration. 
Note that, if the recombination term in Equation 
(48) is omitted (as in Furlanetto et al. 2004), and 
¢ is assumed to be a constant, then every atom in 
the Universe eventually becomes ionized — feo is 
an increasing function of time, and thus eventually 
feou exceeds 1/¢ at every location and reionization 
ends. In reality recombinations are not negligible 
at high enough densities, so as feo increases so 
does the second term in Equation (48). Eventu- 
ally recombinations win, and a few percent of the 
gas mass remains neutral even after reionization, 
locked in the LLSs. 

Figure 11 illustrates this model in practice, fol- 
lowing the original analysis of Kaurov (2016). The 
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Fig. 12 Density map shown as intensity (the values of 
overdensity in the color bar range from 10~? to 100) around 
a massive galaxy in one of the CROC simulations at z = 
9. The red and blue colors show ionized and neutral gas, 
respectively. Filaments remain neutral long after nearby 
voids become fully ionized. 


left panel shows the field of reionization redshifts 
from a numerical simulation of the “Cosmic Reion- 
ization On Computers” (CROC) suite (described 
below). For each pixel, we calculated the moment 
it reached a 90% ionization level weighted by mass 
— the vast majority of the cells cross this threshold 
only once. The top right panel depicts a sample 
of trajectories (or random walks) at z = 8 (about 
the midpoint of reionization) from the Gaussian 
random field initial conditions for this simulation: 
those corresponding to the locations of gas that 
remains neutral (zppy < 8) are colored in blue, 
while the trajectories of cells that become ion- 
ized are shown as red lines. The PDF for neutral 
gas trajectories at a given value of the smoothing 
scale R (using a sharp-k-space filter) is plotted in 
the right bottom panel. The distribution is broad 
and exhibits no sharp features, implying that the 
barrier is “fuzzy”. The ratio between the PDF of 
neutral trajectories and the PDF of all trajecto- 
ries (i.e. the PDF of all overdensities 6) is the 
cumulative barrier — if the barrier was infinitely 
sharp, that ratio would be a Heaviside function 
6[6,(R)—6]. Hence, choosing a value of 6 where this 
ratio is 1/2 serves as a proxy for the barrier. Thus 
defined barrier is shown in the top right panel as 
the green line. 

The fact that the barrier is “fuzzy” should 
not come as a surprise. Figure 12 shows the ion- 
ized bubble around a massive galaxy from one of 
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the CROC simulations at z = 9. Since ionization 
fronts propagate much faster into voids than along 
dense filaments (the I-front speed is proportional 
to ng! =), the spatial distribution of ionized gas 
would not correlate well with the density distri- 
bution on all but the largest spatial scales. Hence, 
using density as a “predictor” for the ionization 
state of the gas is going to be inaccurate. 

Figure 13 shows how the reionization barrier 
evolves towards later epochs in this particular 
realization. By z = 7.3 the barrier reaches the low- 
est gas densities, and ceases to be meaningfully 
defined. At even lower redshifts the distribution 
of neutral gas shifts toward the highest densities, 
which are poorly represented by the correspond- 
ing values in the linear Gaussian field. The right 
panel in the figure shows that the actual non-linear 
PDF of neutral gas has a characteristic overden- 
sity value — A; from Eq. (41), and therefore the 
outside-in model provides an approximate descrip- 
tion for the evolution of the remaining neutral gas. 
It does fail, however, in its main assumption that 
all gas above A; is neutral — in fact, most of the gas 
at any density is ionized, because radiation sources 
are highly clustered, and hence disproportionately 
populate the highest density regions. 


6 Numerical Models of 
Reionization 


6.1 Semi-numerical models 


In this category we place techniques that do not 
require accurate knowledge of the full nonlinear 
matter distribution from numerical simulations. 
Information about the spatial distribution of HII 
regions is obtained directly from the Gaussian ran- 
dom field of the initial conditions or their direct 
derivative, such as the second order perturbation 
theory. 


6.1.1 21cmFAST 


The most widely used in this class of models is per- 
haps the 21cmFAST code (Mesinger et al., 2011). 
In its original formulation it uses the excursion 
set-based formalism of Furlanetto et al. (2004) to 
create a realization of the distribution of ionized 
gas at any given redshift from equations (46) and 
(47). The actual density field used in equation 
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Fig. 13 Left: Same as the bottom right panel of Fig. 11 at four different redshifts. In this particular realization, the lowest 
density gas is ionized at z = 7.3, and at lower redshifts the reionization barrier cannot be meaningfully defined. Right: Same 
PDFs at lower redshifts, but now as function of the total gas density. The distributions of neutral gas have well defined 
peaks (A; from Eq. 41), although they are not cut off sharply at lower densities, so A; is also somewhat “fuzzy”. 


(46) is obtained from the Zel’dovich approxima- 
tion (a historical alias to first order Lagrangian 
perturbation theory). 

One of the most undesirable features of excur- 
sion set-based algorithms is that they do not 
explicitly conserve photons — when two HII bub- 
bles overlap, LyC photons from both bubbles 
are double-counted in setting the ionization field 
(Zahn et al., 2011; Paranjape et al., 2016; Hut- 
ter, 2018). The 21cmFAST code includes a special 
correction to overcome this deficiency (Park et al., 
2021). The latest version of the code is described 
in Murray et al. (2020). It is currently the de- 
facto standard for making large-scale maps of 
various reionization-related fields (shown in Fig. 
14), including creating mock maps for current and 
future 21cm cosmology experiments. 

The excursion set formalism can also be used 
in “reverse”, to model the distribution of neu- 
tral gas — in particular, the distribution of neutral 
islands at the end of reionization. Since we showed 
above that the excursion set formalism loses its 
accuracy close to the overlap, capturing neutral 
islands at this epoch is more involved then just 
“flipping the sign” in the 21cmFAST predictions. 
Such an extension for 21cmFAST has been imple- 
mented in the “IslandFAST” code (Xu et al., 
2017), which agrees well with the results of fully 
coupled simulations. 
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6.1.2 ARTIST-like schemes 


The lack of explicit photon conservation in the 
excursion set-based models has motivated the 
development of alternative approaches that explic- 
itly enforce such conservation. One example is the 
“Asymmetric Radiative Transfer In Shells Tech- 
nique” (ARTIST, Molaro et al., 2019), which 
solves the actual radiative transfer equation by 
ray-tracing (see Section 3.1.1) on a density field 
obtained with the excursion set formalism. To 
make the calculation computationally feasible, 
ARTIST propagates the radiation emitted by each 
source in only a small number of directions, and 
explicitly achieves photon conservation in finite- 
size segments of the solid angle. 

The comparison between ARTIST and excur- 
sion set-based models with the same radiation 
sources shows significant differences in the ion- 
ization history or in the 21cm power spectrum. 
This unfortunately implies that these approaches 
are not robust, and contain internal unidentified 
biases (Molaro et al., 2019). On the other hand, 
one can adjust source parameters such as fesc, 
treating them as fudge factors rather than physi- 
cal quantities, in order to match more closely the 
results of the different algorithms. 

The “Semi-numerical Code for ReionIzation 
with PhoTon Conservation” (SCRIPT, Maity & 
Choudhury, 2022) is an intermediate approach 
between what we call a “semi-numeric technique” 
and a “decoupled simulation”. It can use either 
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Fig. 14 An example of outputs from the latest version of the 21cmFAST code. Maps of several physical quantities 
(such as gas overdensity, neutral fraction, number of recombinations, etc) can be generated along the lightcone. 
(Used by permission from A. Mesinger, ). 


the density field from a dark matter-only (DMO) divided into cells and cells that contain sources 
simulation or the one generated from Gaussian are ionized first, followed by their neighbors, then 
random initial conditions with, say, perturba- neighbors-of-neighbors, etc. We list SCRIPT in 
tion theory. The ionization state is modeled with 
a quasi-ray-tracing algorithm, in which space is 
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Fig. 15 Visualization of the reionization-redshift fields from a full radiative transfer simulation, AMBER, 21cmFAST, and 
“The Reionization on Large Scales” (RLS) scheme from Battaglia et al. (2013). AMBER offers the best approximation to 
the full numerical simulation out of these three approximate techniques. In the units on the x, y axes, the symbol h denotes 
the Hubble constant in units of 100 km s~! Mpc~!. (Used by permission from Hy Trac.) 


this section because it is conceptually very simi- 
lar to ARTIST when used with the density field 
generated from the perturbation theory. 


6.1.3 AMBER 


The “Abundance Matching Box for the Epoch of 
Reionization” (AMBER, Trac et al., 2022) is an 
entirely different approach for modeling cosmic 
dawn. It is based on an “abundance matching” 
hypothesis that there exists a monotonic relation 
between a base field 6(x) and the reionization- 
redshift field at a given location z,¢i(x) (see, e.g., 
example, the left panel of Fig. 11). The functional 
form of Z;e;(8) can be deduced from the assumed 
ionization history. 

The accuracy of such an ansatz depends crit- 
ically on the choice for the base field 6(x). Trac 
et al. (2022) showed that using density as the base 
field is a poor choice, a conclusion that is gener- 
ally consistent with our discussion in Section 5.3. 
A much better choice for 6 is an approximation 
to the ionizing radiation field. The formal solution 
to the radiative transfer equation for the radiation 
angle-averaged intensity J, is 


where 9, is the source function and the 1/(47r?) 
term is the inverse-square law for the flux. Since 
this equation is not a convolution, it cannot be 
easily evaluated on an arbitrary mesh or particle 
distribution and would require using some of the 
methods for solving the radiative transfer equation 
discussed in Section 3.1. 
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There is, however, an ansatz that permits a 
fast — O(N) or O(N log N) — evaluation for an 
arbitrary set of N target points x, and which, on 
average, gives the correct value for the radiation 
field: 


Ty (X,X1) > |x — xi]. 


1 
Amfp (v) 
Equation (49) then becomes a convolution, 


le >e4 Il 


= 3 Si(x1) mip 7) 
oe fe £1 es fe) (50) 


J, (x) 


which can be solved in O(N) operations for N 
spatial locations. 

Figure 15 replicates Fig. 17 from the AMBER 
method paper (Trac et al., 2022), and shows 
a comparison between a full radiative transfer 
simulation, AMBER, 21cmFAST, and an older 
approximate scheme called “The Reionization on 
Large Scales” (RLS, Battaglia et al., 2013). Out 
of the three schemes, AMBER appears to perform 
the best, with 2lcmFAST missing some small scale 
features but otherwise offering a similarly good 
approximation on large scales. 


6.2 DMO+SAM 


In this category we place techniques that rely 
on a dark-matter-only (DMO) simulation com- 
bined with a semi-analytical model (SAM). Such 
an approach is widely used in modeling galaxy for- 
mation, and our focus here is on the small subset 
of such schemes that are specifically designed for 
modeling reionization and, therefore, include an 
algorithm for computing the spatial distribution 
of ionized hydrogen. In the two models discussed 


Fig. 16 Figure 4 from Hutter et al. (2021) showing the 
evolution of the volume-averaged neutral hydrogen frac- 
tion from several ASTRAEUS models and their comparison 
with observational constraints (grey points) — see Hutter 
et al. (2021) for details on the models and the data. The 
good agreement between the different models illustrates 
the robustness of the method, but the discrepancy with the 
observations at z < 5.8 suggests that ASTRAEUS does not 
correctly capture the evolution of the mean free path at 
these epochs. (Used by permission from P. Dayal.) 


below this is done with an excursion set-based 
method. 


6.2.1 ASTRAEUS 


ASTRAEUS (“seminumerical rAdiative tranSfer 
coupling of galaxy formaTion and Reionization 
in N-body dArk mattEr simUlationS”, Hutter 
et al., 2021) is a semi-analytical model (“DEL- 
PHI”, Dayal et al., 2014) built on the “Very 
Small MultiDark Planck” DMO simulation from 
the MultiDark simulation suite (Klypin et al., 
2016). The key parameters of the simulation are 
listed in Table 1. The ionization state of the IGM 
is modeled by an independent excursion set code 
developed by Hutter (2018). With a volume size 
of © 230 Mpc, ASTRAEUS simulations are large 
enough to produce a converged distribution of 
ionized bubbles during reionization (Iliev et al., 
2014). 

Figure 16 shows an example from ASTRAEUS 
modeling. The global reionization history in 
ASTRAEUS depends only weakly on the specifics 
of the adopted radiative feedback model, in agree- 
ment with many other studies. On a negative side, 
ASTRAEUS fails to model the ionization state of 
the IGM in the post-reionization era at z < 6, 
and thus cannot yet be used for making synthetic 
quasar spectra for studying transmission spikes, 
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Fig. 17 Figure 4 from Mutch et al. (2016) showing the 
electron scattering optical depth to reionization, Te, as a 
function of redshift for several DRAGONS models - see 
Mutch et al. (2016) for details on the simulations and 
observational constraints. (Used by permission from S. 
Mutch.) 


dark gaps, and the Lyman-a forest. On a positive 
side, the large volumes simulated by ASTRAEUS 
allow one to explore sample variance in current 
and future observational surveys. For example, 
Ucci et al. (2021) found large variations in the 
galaxy luminosity functions at Myy < —17 that 
will soon be measured by the JWST Advanced 
Extragalactic Survey (JADES). 


6.2.2 DRAGONS 


DRAGONS (“Dark-ages Reionization And 
Galaxy Observables from Numerical Simula- 
tions”, Poole et al., 2016; Angel et al., 2016; 
Mutch et al., 2016) adopt a similar scheme 
to ASTRAEUS. It uses the DMO simulation 
“Tiamat” that was specifically designed for the 
DRAGONS project. While 12 times smaller in 
volume and about 6 times smaller in particle 
number than the “Very Small MultiDark Planck” 
simulation used in ASTRAEUS, Tiamat has more 
than 10 times better spatial resolution. 

The power of semi-numerical schemes like 
ASTRAEUS and DRAGONS is in their flexibility. 
After the backbone DMO simulation is completed, 
running the semi-analytical prescriptions on the 
halo catalogs is straightforward and fast, and can 
therefore be used for an efficient exploration of a 
large parameter space. As an example, we show in 
Fig. 17 several DRAGONS models with different 
stellar feedback recipes (not to be confused with 
the radiative feedback models shown in Fig. 16 


for ASTRAEUS). Models without the stellar feed- 
back predict too early reionization, but, of course, 
the feedback is well established to be the key 
in modeling galaxy formation, so models with- 
out feedback are of only academic interest. With 
strong enough feedback, star formation in galax- 
ies “self-regulates” (Dobbs et al., 2011; Agertz 
et al., 2013; Hopkins et al., 2013; Agertz & 
Kravtsov, 2015; Benincasa et al., 2016; Hopkins 
et al., 2018; Orr et al., 2018; Semenov et al., 
2017; Faucher-Giguére, 2018; Grudi¢ et al., 2018; 
Semenov et al., 2018; Li et al., 2018; Semenov 
et al., 2019; Meng et al., 2019; Barrera-Ballesteros 
et al., 2021; Grudié et al., 2021), and the main 
parameter controlling the timing of reionization is 
the escape fraction. 


6.2.3 1D radiative transfer 


An interesting simplified approach to implement- 
ing radiative transfer on top of a DMO simulation 
was pioneered by Thomas & Zaroubi (2008) and 
later developed into the GRIZZLY (not an abbre- 
viation) code (Ghara et al., 2015, 2018). Instead 
of utilizing the full ray-tracing, only one ray is 
used for each source, with the density profile along 
the ray being the spherical average over all direc- 
tions. The radiative transfer is solved along that 
ray in “1D”, so all ionized regions around sources 
are spherical. This approach has been shown to 
be a useful approximation for modeling the 21 cm 
emission during reionization (Ghara et al., 2018). 
It should be clear from our discussion in Section 5, 
however, that such approach fails to capture the 
“outside-in” stage of reionization. 


6.3 Partially coupled simulations 


Cosmological simulations vary widely in their 
setups, input physics, and numerical approaches. 
In this review we deliberately limit ourselves to 
simulations that model reionization, i.e. include 
cosmological radiative transfer as one of their 
main physics packages. A number of recent com- 
putational projects have focused on modeling high 
redshift galaxy formation without accounting for 
the thermal and ionization history of their envi- 
ronment, such as FLARES (Lovell et al., 2021; 
Vijayan et al., 2021, 2022), FIRE (Ma et al., 
2020b; Liang et al., 2021), or FirstLight (Ceverino 
et al., 2017). Because they do not address the 
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reionization process itself, they are somewhat out- 
side the scope of the present review. An excellent 
recent review of galaxy formation simulations is 
given by Vogelsberger et al. (2020). 

An intermediate stage between a “semi- 
numerical” technique like DMO+SAM and a fully 
coupled cosmological hydrodynamic simulation is 
one where not all the physics components of a sim- 
ulation are actually fully coupled. The classical 
example of such an approach are C2-Ray simula- 
tions (Iliev et al., 2006b; Mellema et al., 2006), 
which include full radiative transfer but assume 
that the hydrogen density tracks the dark mat- 
ter field of a large DMO run. In this method, 
the gas dynamics does not respond to heating 
by UV radiation, and physical effects such as the 
suppression of gas accretion and condensation in 
sufficiently low-mass halos cannot be captured 
directly, although they may be modeled with addi- 
tional, approximate schemes. These simulations 
also do not resolve the scales that are relevant for 
star formation, and model radiation sources with 
a semi-analytic approach. 

A different example of partially coupled sim- 
ulations can be found in Bauer et al. (2015), 
where two different radiative transfer solvers are 
run in post-processing on outputs from the Illus- 
tris simulation (Vogelsberger et al., 2014). Illustris 
aims at studying the processes of galaxy for- 
mation and evolution in the Universe with a 
comprehensive physical model that includes feed- 
back from supernova explosions and accreting 
supermassive black holes, and radiative transfer 
in post-processing complements the original sim- 
ulation with maps of ionized/neutral IGM gas 
throughout the reionization process. The limita- 
tion of this approach is that the gas dynamics does 
not respond to the spatially-varying and time- 
dependent radiation field. An analogous approach 
using the MassiveBlack-II simulation (Khandai 
et al., 2015a) and the radiative transfer code 
CRASH was presented by Eide et al. (2018, 2020); 
Ma et al. (2022). Here, the key advantage is the 
adopted frequency range, which is wider than 
modeled by most other approaches and extend to 
soft X-rays and the associated secondary ioniza- 
tions. The CRASH+MassiveBlack-II set of simu- 
lation also systematically explored the impacts of 
multiple types of radiation sources — in addition 
to normal stars and quasars, they also considered 
binary stars, shock-heated ISM, and X-ray bursts. 


The radiative transfer solver has also been used in 
post-processing for the modeling, for example, of 
the escape fraction from individual halos in cosmo- 
logical simulations (c.f. Ma et al., 2020b; Kostyuk 
et al., 2022, for recent studies). As we mentioned 
at the beginning of this subsection, such models 
are beyond the scope of this review. 

Partially coupled simulations with a primary 
focus on galaxy formation at high redshift and 
not on the process of reionization per se may also 
include an approximate radiative transfer algo- 
rithm rather than an actual numerical solver. The 
poster child examples of such schemes are Blue- 
Tides (Huang et al., 2018) and Astrid (Bird et al., 
2022). These simulations combine the physics 
package of the MassiveBlack-II run (Khandai 
et al., 2015b) with the “The Reionization on Large 
Scales” approach from Battaglia et al. (2013) that 
served as a precursor to the AMBER scheme 
already described. They are able to account, albeit 
approximately, for the ionization history of the 
galaxy environment, while offering a huge boost 
to computational efficiency. 

Yet another flavor of a partially decoupled 
scheme can be envisioned as an efficient approach 
for a multi-fold extension of the simulation size. 
Most of the methods described here as well as the 
fully coupled simulations discussed below resolve 
the low density IGM on scales of order 100 
comoving kpc or better. Such resolution, com- 
parable to the Doppler smoothing scale of gas 
at 104K, is the minimum required to capture 
the IGM density fluctuations that give origin 
to the Lyman-a forest in the spectra of distant 
quasars. A uniform grid simulation with, say, 
50 kpc resolution in a 500 Mpc box would require 
10,000% cells, a size that is currently achievable 
on modern 100 petaflops-level platforms with, 
e.g., the GPU-native cosmological hydrodynam- 
ics code Cholla (Schneider & Robertson, 2015). 
At 50kpc resolution, however, no plausible model 
for the physics of galaxy formation can be con- 
structed. This is similar, for example, to C2-Ray 
simulations that have to rely on semi-analytical 
schemes for including radiation sources. A “two- 
tier” simulation approach may be desirable, where 
a high-resolution small-box simulation is used to 
inform a coarse-resolution large-box one. Large 
volume simulations — which cannot track the inte- 
rior structure of dark matter halos — may then 
implement the physics of galaxy formation with 
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Fig. 18 The EBM model for the star formation rate 
(SFR) as a function of virial mass Myir. The upper left 
panel shows the two-dimensional distribution of SFRs with 
Myjiy for all galaxies in a subset of CROC simulations, with 
the color scale showing the number of simulated galaxies 
at each [SFR, Myir] location. The lower left panel shows 
the EBM model prediction I'(SFR; 6’) for the distribution 
of SFR with M,yi; given the EBM internal parameters 6’. 
The upper and lower right panels show the residuals and 
the outliers between the simulated CROC galaxies and the 
EBM model predictions. These outliers represent less than 
5% of the simulated CROC galaxies. (Used by permission 
from Ryan Hausen.) 


a approximate model that recovers the mean 
trends for galaxy baryonic properties predicted by 
more detailed calculations. A recent example of 
this technique is offered by Hausen et al. (2022), 
who trained the Explainable Boosting Machines 
(EBM) machine learning algorithm to assign stel- 
lar masses and star formation rates (SFRs) to the 
host dark matter halos based on nearly 6 million 
galaxies simulated by the fully coupled “Cosmic 
Reionization On Computers” (CROC) project. 
Figure 18 shows the main result from that work: a 
comparison of the simulated SFRs in the original 
CROC simulations versus the values predicted by 
machine learning as a function of halo virial mass 
M,ir. The EBM model is highly predictive, fail- 
ing to capture only the most extreme outliers in 
SFR at a fixed M,;,. Through this approach, the 
physics of baryonic galaxy formation can be con- 
nected to the properties of dark matter halos and 
implemented as a “sub-grid” prescription in cos- 
mological hydrodynamics simulations that do not 


resolve the small scale details of star formation 
and feedback, while at the same time capturing 
the variations of SFR in halos of the same mass 
due to environmental effects and different prior 
accretion histories. 


6.4 Fully coupled simulations 


Self-consistent, fully coupled simulations, often 
considered the ultimate theoretical model for a 
given process, can only be trusted as long as 
they include all the relevant physics (such as 
gravity, gasdynamics, star formation, stellar and 
AGN feedback, radiation transport) with suffi- 
cient precision and maintain numerical effects 
under control. 

Table 1 summarizes some of the most recent 
fully coupled simulations (together with several of 
the approximate methods discussed above). Cos- 
micDawn (CoDa) (Ocvirk et al., 2016; Aubert 
et al., 2018; Ocvirk et al., 2020; Lewis et al., 
2020, 2022), CROC, and Thesan are similar in 
their computational volume (~ 100 Mpc). They 
differ somewhat in spatial and mass resolution, 
but all fall into the class of simulations that do 
not resolve the scale heights of galactic disks (and 
therefore model galaxies in “2D”), and can be 
directly compared to each other. SPHINX simu- 
lations (Rosdahl et al., 2018; Katz et al., 2019, 
2021), by contrast, focus on resolving the actual 
vertical structure of star-forming galaxies, achiev- 
ing much higher spatial resolution at the expense 
of being unable to model the global reionization 
history (a 20 comoving Mpc box contains only 4 
L, galaxies on average). 

The target goal of reaching a ~ 100 Mpc sim- 
ulation box is dictated by the desire to replicate 
a representative region of the Universe. At z ~ 
7, the correlation length of galaxies is around 
10 Mpc (Barone-Nugent et al., 2014) and is weakly 
luminosity dependent. The number density of L, 
galaxies is also about 1 per 10 Mpc. Hence, a vol- 
ume of ~ 100 Mpc on a side contains around 1000 
L, galaxies and has the rms correlation function at 
half the box size of (50/10)~1° = 0.08. Whether 
~ 100 Mpc simulations actually converge on some 
of the key features of the reionization process, 
such as the size distribution of ionized bubbles, 
is presently unclear. Some earlier C2-Ray simula- 
tions (Iliev et al., 2014) found convergence only 
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in = 250Mpc boxes, while in CROC simula- 
tions (which explicitly match at z < 6 the LyC 
mean free path determined by the abundance of 
LLSs) the bubble size distribution appears to have 
converged by z = 7 (Gnedin & Kaurov, 2014). 

Figure 19 from Bouwens et al. (2022) shows 
a comparison between the galaxy UV luminos- 
ity function observed at different redshifts and 
several theoretical predictions. The observational 
data come primarily from the Hubble Frontier 
Fields, and are subject to systematic uncertain- 
ties from lensing modeling (Bouwens et al., 2017, 
2022) that rapidly increase for galaxies fainter 
than Myy > —15. At brighter magnitudes most 
theoretical models match the data well,’ but they 
differ widely at Myy = —14 mag, a faint lumi- 
nosity regime that will soon be probed by JWST 
observations. Theorists are eagerly waiting for 
JWST deep field data as few, if any, of the cur- 
rent models will survive these ground-breaking 
sensitive measurements. 

Since the goal of fully coupled simulations is to 
actually model the reionization process, their per- 
formance must be measured against that metric. 
In Figure 20 we show the most basic prediction 
for any simulation of the epoch of reionization, the 
evolution of the globally volume-averaged neutral 
hydrogen fraction.® While the latest fully coupled 
simulations do produce neutral fractions similar 
to the observed values, they do not do so at the 
level of precision required by the observations. 
The only simulation project that does match the 
observations accurately is CROC, but this is by 


construction, as the data measurements were actu- 


ally used to calibrate the effective escape fraction 
from the radiation sources. The unsatisfactory 
level of agreement between data and theory points 
to the direction where efforts in designing the 
next generation of cosmological simulations should 
be aimed at — i.e. on significantly improving the 
modeling of intergalactic gas. A number of recent 
observations have provided measurements of vari- 
ous properties of the post-reionization IGM, such 
as the distribution of mean Lyman-a opacities 


“Notice, however, that CROC results are shown for earlier, 
smaller boxes. CROC underpredicts luminosities and stellar 
masses of super-L, galaxies in its largest, 117 Mpc boxes (Zhu 
et al., 2020). 

8 Predicting the mass-weighted neutral hydrogen fraction is 
much harder, as after reionization this is dominated by the 
residual neutral component locked in the Damped Lyman-a 
Systems. 


Table 1 Simulation-based models and simulations of reionization. 


Suite Name Nobox” | Npart? | Box size © | Resolution? | Code ° 
SAM+DMO 
ASTRAEUS 1 3840°% 230 Mpc 530 pe ART (DMO) 
DRAGONS 1 2160 100 Mpc 50 pe GADGET-2 (DMO) 
Partially Coupled Simulations 
Astrid 1 5500% 370 Mpc 115 pe GADGET-3 (SPH) + semi-analytical RT 
MassiveBlack-II 1 17923 | 142 cMpc 440 pe GADGET-3 (SPH) + RT in post-processing 
C2-Ray 1 | 69123 | 714Mpe | 6.3kpe | PM (DMO) + C2-Ray (RT) 
Illustris 1 18203 | 107 Mpc 72 pe AREPO (MM) + RT in post-processing 
Fully Coupled Simulations 
CoDa 1 81923 94 Mpc 1.7 kpe RAMSES-CUDATON (uniform grid) 
CROC 6 | 20483 | 117 Mpe 100 pe | ART (AMR) 
SPHINX 1 1024°% 20 Mpc 10 pe RAMSES-RT (AMR) 
Thesan 1 21002 96 Mpc 300 pc AREPO-RT (MM) 
“We list only the number of the largest box runs. Several projects also completed a number of simulations with 8 times smaller 


Npart- 


°In comoving units. 


’Equivalent number of dark matter particles in the largest simulation. 


4Spatial resolution is given in proper units as an effective grid code cell size; for particle codes the conversion between the 
gravitational softening and the effective cell size is given in Mansfield & Avestruz (2021). For simulations that maintained their 


resolution in comoving units, the resolution is quoted at z = 6. 


°“AMR=Adaptive Mesh Refinement; SPH=Smooth Particle Hydrodynamics; MM=Moving Mesh; RT=Radiative Transfer. 


along skewers of fixed length (Becker et al., 2015; 
Bosman et al., 2018; Eilers et al., 2018; Yang 
et al., 2020), the distribution of “dark gaps” — 
continuous spectral regions in distant quasar spec- 
tra where the trasmitted flux is below a specified 
threshold (Zhu et al., 2020, 2021), and the cross- 
correlation between quasar absorption spectra and 
properties of galaxies along the same sightline 
(Meyer et al., 2019, 2020). None of the existing 
fully coupled simulations are able to match all of 
these observational constraints. 


7 Concluding Remarks 


The epoch of cosmic reionization is a vast science 
frontier discovery area where detailed modeling 
and complex computational tools are required 
to support the anticipated observational break- 
throughs. Since the very first paper on reion- 
ization by Arons & McCray (1970), this field 
has remained heavily theory dominated until the 
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very recent times. The prospects for improving 
our understanding of this era over the next few 
years are outstanding, and a number of mas- 
sive datasets from new ground- and space-based 
instruments and facilities, like the Dark Energy 
Spectroscopic Instrument (DESI), the Subaru 
Hyper Suprime-Cam (HSC), the Extremely Large 
Telescope (ELT), the Atacama Large Millimeter 
Array (ALMA), the Hydrogen Epoch of Reioniza- 
tion Array (HERA), the Square Kilometre Array 
(SKA), JWST and Euclid, are poised to revolu- 
tionize our understanding of the dawn of galaxies, 
the cosmic ionizing photon budget, the physics of 
the early IGM and galaxy-IGM interactions, and 
promise an era of precision reionization studies. 
It is generally expected that, by perform- 
ing throughout analyses that comprise the latest 
observational results, elaborate theoretical mod- 
eling, state-of-the-art simulations that incorpo- 
rate realistic baryonic physics over a large range 
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Fig. 19 Comparison (Bouwens et al., 2022) of the observed galaxy UV luminosity functions (gray bands) with model predictions 
from DRAGONS, CROC, “Technicolor Dawn” (F17, Finlator et al., 2018), 2lcmFast (P19, Park et al., 2019), a semi-analytical 
model (Y19, Yung et al., 2019), “Cosmic Dawn” (020, Ocvirk et al., 2020), ASTRAEUS, Delphi (an earlier version of ASTRAEUS), 
and FirstLight (Ceverino et al., 2017). (Used by permission from R. Bouwens.) 


of spatial scales, and statistical inference, astro- 
physicists will be able over the next decade to: 
1) unravel the nature of the first astrophysical 
sources of radiation and heating, and investi- 
gate how they interacted with their environment; 
2) trace the overall progression of the cosmic 
reionization process and establish its impact on 
observations of the early Universe — from the large- 
scale signal of neutral hydrogen in the redshifted 
21cm line at z = 10 to the fluctuating Lyman-a 
opacity of the z ~ 5.5 IGM; 3) reveal the ther- 
mal history of cosmic gas through the imprint left 
by hydrogen and helium reionization heating on 
the power spectrum of the Lyman-a forest and 
other observables; and 4) model the range of cos- 
mological physics that influence the small-scale 
properties of the cosmic web. 

The precision of future observational data will 
exceed the precision of current theoretical models 
(this is already happening, as can be seen in Fig. 
20). Theorists must therefore face the challenge 
of significantly increasing the physical fidelity and 
numerical precision of their models, and enhance 
their predictive power. The world of theoretical 
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modeling of reionization is very diverse, as we 
strove to demonstrate above. Such a diversity is 
useful during the formative years of any field, but 
as the field matures one expects many of the sim- 
pler, less accurate methods and tools to become 
obsolete. This is very likely going to be the case 
in reionization studies as well. It is our belief then 
that this “living” review will lead a very active 
life indeed, with many of the schemes described 
above falling out of use. Hopefully, this article will 
be of help to the next generation of theorists and 
computational scientists, charting the directions 
where a beginning graduate student or a post- 
doc switching fields may make a large and visible 
impact. 
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